diff --git a/.flake8 b/.flake8 index 969b66bdb..696d05318 100644 --- a/.flake8 +++ b/.flake8 @@ -1,7 +1,9 @@ [flake8] -ignore = E402 C901 E501 E265 +ignore = E402 exclude = .eggs, *.egg,build, src/mrsimulator/__init__.py filename = *.pyx, *py max-line-length = 88 -max-complexity = 10 +max-complexity = 12 select = C,E,F,W,N8 +count = True +statistics = True diff --git a/conftest.py b/conftest.py index af04953be..5a81ddcfd 100644 --- a/conftest.py +++ b/conftest.py @@ -3,6 +3,8 @@ import matplotlib import matplotlib.pyplot as plt +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo import numpy as np import pytest from mrsimulator import Simulator @@ -28,6 +30,8 @@ def add_site(doctest_namespace): doctest_namespace["st"] = SymmetricTensor doctest_namespace["pprint"] = pprint doctest_namespace["Isotope"] = Isotope + doctest_namespace["sp"] = sp + doctest_namespace["apo"] = apo site1 = Site( isotope="13C", diff --git a/docs/api_py/fitting.rst b/docs/api_py/fitting.rst index a778e5e33..38ca18bd9 100644 --- a/docs/api_py/fitting.rst +++ b/docs/api_py/fitting.rst @@ -5,5 +5,5 @@ Fitting Utility Functions .. currentmodule:: mrsimulator.spectral_fitting -.. autofunction:: make_fitting_parameters -.. autofunction:: min_function +.. autofunction:: make_LMFIT_parameters +.. autofunction:: LMFIT_min_function diff --git a/docs/api_py/operations/operation_documentation.rst b/docs/api_py/operations/operation_documentation.rst new file mode 100644 index 000000000..1049b437f --- /dev/null +++ b/docs/api_py/operations/operation_documentation.rst @@ -0,0 +1,15 @@ + + +Documentation +------------- + +.. currentmodule:: mrsimulator.signal_processing + +.. autoclass:: Scale +.. autoclass:: IFFT +.. autoclass:: FFT + +.. currentmodule:: mrsimulator.signal_processing.apodization + +.. autoclass:: Gaussian +.. autoclass:: Exponential diff --git a/docs/api_py/operations/operations_summary.rst b/docs/api_py/operations/operations_summary.rst new file mode 100644 index 000000000..e075071b1 --- /dev/null +++ b/docs/api_py/operations/operations_summary.rst @@ -0,0 +1,49 @@ +.. _operations_api: + +Operations +========== + +Generic operations +------------------ + +.. currentmodule:: mrsimulator.signal_processing + +Import the module as + +.. doctest:: + + >>> import mrsimulator.signal_processing as sp + +.. rubric:: Operation Summary + +The following list of operations apply to **all dependent variables** within the +simulation data object. + +.. autosummary:: + :nosignatures: + + ~Scale + ~IFFT + ~FFT + +Apodization +----------- + +.. currentmodule:: mrsimulator.signal_processing.apodization + +Import the module as + +.. doctest:: + + >>> import mrsimulator.signal_processing.apodization as apo + +.. rubric:: Operation Summary + +The following list of operations apply to **selected dependent variables** within +the simulation data object. + +.. autosummary:: + :nosignatures: + + ~Gaussian + ~Exponential diff --git a/docs/api_py/py_api.rst b/docs/api_py/py_api.rst index e382ebe7a..cd76b8d43 100644 --- a/docs/api_py/py_api.rst +++ b/docs/api_py/py_api.rst @@ -15,6 +15,8 @@ Python-API References method methods fitting + signal_processing + operations/operations_summary .. parameterized_tensor .. interaction diff --git a/docs/api_py/signal_processing.rst b/docs/api_py/signal_processing.rst new file mode 100644 index 000000000..9cfcdf9f7 --- /dev/null +++ b/docs/api_py/signal_processing.rst @@ -0,0 +1,15 @@ +.. _signal_processing_api: + +Signal Processing +================= + +.. currentmodule:: mrsimulator.signal_processing + +.. autoclass:: SignalProcessor + :show-inheritance: + + .. rubric:: Method Documentation + + .. automethod:: parse_dict_with_units + .. automethod:: to_dict_with_units + .. automethod:: apply_operations diff --git a/docs/api_py/simulator.rst b/docs/api_py/simulator.rst index 0e4c6f1ea..af27034e4 100644 --- a/docs/api_py/simulator.rst +++ b/docs/api_py/simulator.rst @@ -20,4 +20,3 @@ Simulator .. automethod:: run .. automethod:: save .. automethod:: load - .. automethod:: apodize diff --git a/docs/index.rst b/docs/index.rst index ec1092fde..859a9353a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -152,6 +152,7 @@ Our current objectives for the future are the following getting_started using_mrsimulator_objects configuring_simulator + signal_processing mrsim_IO benchmark auto_examples/index diff --git a/docs/requirements.rst b/docs/requirements.rst index a23ceec34..709f967c1 100644 --- a/docs/requirements.rst +++ b/docs/requirements.rst @@ -15,7 +15,7 @@ Package dependencies - `requests>=2.22 `_ - cython>=0.29.14 - `matplotlib>=3.1 `_ for figures and visualization, -- `csdmpy>=0.3 `_ +- `csdmpy>=0.3.1 `_ - `pydantic>=1.0 `_ - monty>=2.0.4 diff --git a/docs/signal_processing.rst b/docs/signal_processing.rst new file mode 100644 index 000000000..80528146f --- /dev/null +++ b/docs/signal_processing.rst @@ -0,0 +1,331 @@ + +Post Simulation Signal Processing +================================= +.. sectionauthor:: Maxwell C. Venetos + + +After running a simulation, you may often find yourself in need of some post-simulation +signal processing operations. For example, you may want to scale the intensities to +match the experiment, add line-broadening, or simulate signal artifact such as sinc +wiggles. There are many signal-processing libraries, such as Numpy and Scipy, that you +may use to accomplish this. Although, in NMR, certain operations like applying +line-broadening, is so regularly used that it soon becomes inconvenient to having to +write your own set of code. For this reason, the ``mrsimulator`` package offers some +frequently used signal-processing operations. + +The following section will demonstrate the use of the +:class:`~mrsimulator.signal_processing.SignalProcessor` class in applying various +operations to the simulation data. + +**Setup for the figures** + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> import matplotlib as mpl + >>> import matplotlib.pyplot as plt + ... + >>> # global plot configuration + >>> font = {"weight": "light", "size": 9} + >>> mpl.rc("font", **font) + >>> mpl.rcParams["figure.figsize"] = [4.5, 3] + + +Simulating spectrum +------------------- +Please refer to the :ref:`using_objects` for a detailed description +of how to set up a simulation. Here, we will create a hypothetical simulation from two +single-site spin-systems to illustrate the use of the post-simulation signal processing +module. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> from mrsimulator import Simulator, SpinSystem, Site + >>> from mrsimulator.methods import BlochDecaySpectrum + ... + >>> # **Step 1** Create the site and spin system objects + >>> site1 = Site( + ... isotope="29Si", + ... isotropic_chemical_shift=-75.0, # in ppm + ... shielding_symmetric={"zeta": -60, "eta": 0.6}, # zeta in ppm + ... ) + >>> site2 = Site( + ... isotope="29Si", + ... isotropic_chemical_shift=-80.0, # in ppm + ... shielding_symmetric={"zeta": -70, "eta": 0.5}, # zeta in ppm + ... ) + ... + >>> sites = [site1, site2] + >>> labels = ["Sys-1", "Sys-2"] + >>> spin_systems = [SpinSystem(name=l, sites=[s]) for l, s in zip(labels, sites)] + ... + >>> # **Step 2** Create a Bloch decay spectrum method. + >>> method_object = BlochDecaySpectrum( + ... channels=["29Si"], + ... magnetic_flux_density=7.1, # in T + ... rotor_angle=54.735 * 3.1415 / 180, # in rads + ... rotor_frequency=780, # in Hz + ... spectral_dimensions=[ + ... { + ... "count": 2048, + ... "spectral_width": 25000, # in Hz + ... "reference_offset": -5070, # in Hz + ... "label": "$^{29}$Si resonances", + ... } + ... ], + ... ) + ... + >>> # **Step 3** Create the simulation and add the spin system and method objects. + >>> sim = Simulator() + >>> sim.spin_systems = spin_systems + >>> sim.methods = [method_object] + ... + >>> # **Step 4** Simulate the spectra. + >>> sim.run() + +The plot the spectrum is shown below. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> ax = plt.subplot(projection="csdm") # doctest: +SKIP + >>> ax.plot(sim.methods[0].simulation, color="black", linewidth=1) # doctest: +SKIP + >>> ax.invert_xaxis() # doctest: +SKIP + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP + +Post-simulating processing +-------------------------- + +Signal processing is a series of operations that are applied to the dataset. In this +workflow, the result from the previous operation becomes the input for the next +operation. In the ``mrsimulator`` library, we define this series as a list of operations. + +Setting a list of operations +'''''''''''''''''''''''''''' + +All signal processing operations are located in the `signal_processing` module of the +``mrsimulator`` library. Within the module is the `apodization` sub-module. An +apodization is a point-wise multiplication operation of the input signal with the +apodizing vector. Please read our :ref:`operations_api` documentation for a complete +list of operations. + +Import the module and sub-module as + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> import mrsimulator.signal_processing as sp + >>> import mrsimulator.signal_processing.apodization as apo + +In the following example, we show the application of a single operation—-convoluting +the frequency spectrum with a Gaussian lineshape, that is, simulating a Gaussian +line-broadening--using the :class:`~mrsimulator.signal_processing.SignalProcessor` +class. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> # list of processing operations + >>> post_sim = sp.SignalProcessor( + ... operations=[ + ... sp.IFFT(), apo.Gaussian(sigma=100), sp.FFT() + ... ] + ... ) + +The required attribute of the ``SignalProcessor`` class, `operations`, holds the list of +operations that gets applied to the input dataset. The above set of operations is for a +frequency domain input signal undergoing a Gaussian convolution of 100 Hz. In this scheme, +the operations list will first perform an inverse Fourier Transform to convert +the frequency domain signal to the time domain. Next, the time domain signal is apodized +by a Gaussian function with a broadening factor of 100 Hz, followed by a forward Fourier +transformation transforming the signal back to the frequency domain. + +.. note:: + For almost all NMR spectrum, the post-simulation processing is a convolution, including + the line-broadening. The convolution theorem states that under suitable conditions, the + Fourier transform of a convolution of two signals is the pointwise product of their + Fourier transforms. + + +Applying operation to the spectrum +'''''''''''''''''''''''''''''''''' + +To apply the above list of operations to the simulation/input data, use the +:meth:`~mrsimulator.signal_processing.SignalProcessor.apply_operations` method of the +``SignalProcessor`` instance as follows, + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) + +The `data` is the required argument of the `apply_operations` method, whose value is a +CSDM object holding the dataset. The variable `processed_data` holds the output, that is, +the processed data. The plot of the processed signal is shown below. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> ax = plt.gca(projection="csdm") # doctest: +SKIP + >>> ax.plot(processed_data, color="black", linewidth=1) # doctest: +SKIP + >>> ax.invert_xaxis() # doctest: +SKIP + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP + + +Applying operation to the sub-spectra +''''''''''''''''''''''''''''''''''''' + +.. spectrum and follow up by decomposing the spectrum and processing each signal +.. independently. +.. The above code resulted in the same processing to be applied +.. to both signals because in the simulation the signals were not +.. seperated. + +It is not uncommon for the NMR spectrum to compose of sub-spectrum, from different +sites/systems, exhibiting differential relaxations, and therefore, have different +extents of line-broadening. The reason for this differential relaxation behavior is +not the focus of this sub-section. Here, we show how one can simulate such spectra +using the operations list. + +Before we can move forward, you will first need to identify these sub-systems and +simulate individual spectra for these systems. In this example, we will treat the two +spin-systems as the two different spin environments exhibiting different +relaxations/line-broadening. To simulate the sub-spectrum from the individual +spin-systems, modify the value of the :attr:`~mrsimulator.Simulator.config` attribute +as follows, and re-run the simulation. +Refer to the :ref:`config_simulator` section for further details. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> sim.config.decompose_spectrum = "spin_system" + >>> sim.run() + +.. Note, in the previous example, both sites/spin-systems got the same extent of Gaussian +.. line-broadening. The following example illustrates how you can apply you might want to apply a different set of +.. In order to apply different processes to each signal, +.. we must set the simulation config to decompose the spectrum. +.. Steps 1-3 will be the same and we will start at step 4. +.. # +.. **Step 4** Decompose spectrum and run simulation. +.. sim.config.decompose_spectrum = "spin_system" +.. sim.run() +.. plt.xlabel("$^{29}$Si frequency / ppm") +.. plt.xlim(x.value.max(), x.value.min()) +.. plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) + +The above code generates two spectra, each corresponding to a spin-system. +The plot of the spectra is shown below. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> ax = plt.gca(projection="csdm") # doctest: +SKIP + >>> ax.plot(sim.methods[0].simulation) # doctest: +SKIP + >>> ax.invert_xaxis() # doctest: +SKIP + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP + +Because the simulation is stored as a CSDM [#f1]_ object, each sub-spectrum is a +dependent-variable of the CSDM object, sharing the same frequency dimension. +When using the list of the operations, you may selectively apply a given operation to a +specific dependent-variable by specifying the index of the corresponding +dependent-variable as an argument to the operation class. Note, the order of the +dependent-variables is the same as the order of the spin-systems. Use the `dep_var_indx` +argument of the operation to specify the index. Consider the following list of +operations. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> post_sim = sp.SignalProcessor( + ... operations=[ + ... sp.IFFT(), # convert to time-domain + ... apo.Gaussian(sigma=50, dep_var_indx=0), + ... apo.Exponential(FWHM=200, dep_var_indx=1), + ... sp.FFT(), # convert to frequency-domain + ... ] + ... ) + +The above operations list first applies an inverse Fourier transformation, +followed by a Gaussian apodization on the dependent variable at index 0 (spin-system +labeled as `sys1`), followed by an Exponential apodization on the dependent +variable at index 1 (spin-system labeled as `sys2`), and finally a forward Fourier +transform. Note, the FFT and IFFT operations apply on all dependent-variables. + +As before, apply the operations with the +:meth:`~mrsimulator.signal_processing.SignalProcessor.apply_operations` method. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) + +The plot of the processed spectrum is shown below. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> ax = plt.gca(projection="csdm") # doctest: +SKIP + >>> ax.plot(processed_data, alpha=0.75) # doctest: +SKIP + >>> ax.invert_xaxis() # doctest: +SKIP + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP + + +Serializing the operations list +------------------------------- + +You may also serialize the operations list using the +:meth:`~mrsimulator.signal_processing.SignalProcessor.to_dict_with_units` +method, as follows + +.. doctest:: + + >>> from pprint import pprint + >>> pprint(post_sim.to_dict_with_units()) + {'operations': [{'dim_indx': 0, 'function': 'IFFT'}, + {'dep_var_indx': 0, + 'dim_indx': 0, + 'function': 'apodization', + 'sigma': '50.0 Hz', + 'type': 'Gaussian'}, + {'FWHM': '200.0 Hz', + 'dep_var_indx': 1, + 'dim_indx': 0, + 'function': 'apodization', + 'type': 'Exponential'}, + {'dim_indx': 0, 'function': 'FFT'}]} + +.. [#f1] Srivastava, D. J., Vosegaard, T., Massiot, D., Grandinetti, P. J., + Core Scientific Dataset Model: A lightweight and portable model and + file format for multi-dimensional scientific data, PLOS ONE, + **15**, 1-38, (2020). + `DOI:10.1371/journal.pone.0225953 `_ diff --git a/docs/using_mrsimulator_objects.rst b/docs/using_mrsimulator_objects.rst index 27177186f..d70bf9f6e 100644 --- a/docs/using_mrsimulator_objects.rst +++ b/docs/using_mrsimulator_objects.rst @@ -2,9 +2,6 @@ .. _using_objects: -.. .. image:: https://mybinder.org/badge_logo.svg -.. :target: https://mybinder.org/v2/gh/DeepanshS/mrsimulator/master?filepath=jupyternotebooks%2F - =================================================== Getting started with ``mrsimulator``: Using objects =================================================== diff --git a/environment-dev.yml b/environment-dev.yml index 38b41fb8b..41b74c6aa 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -13,7 +13,7 @@ dependencies: - pip - pip: - monty>=2.0.4 - - csdmpy>=0.3 + - csdmpy>=0.3.1 - pydantic>=1.0 - typing-extensions>=3.7 - flake8 diff --git a/environment.yml b/environment.yml index f1fbd1221..d9d200381 100644 --- a/environment.yml +++ b/environment.yml @@ -12,6 +12,6 @@ dependencies: - pip - pip: - monty>=2.0.4 - - csdmpy>=0.3 + - csdmpy>=0.3.1 - pydantic>=1.0 - typing-extensions>=3.7 diff --git a/examples/1D_simulation/plot_0_Wollastonite.py b/examples/1D_simulation/plot_0_Wollastonite.py index f5ab6845f..2c7da0348 100644 --- a/examples/1D_simulation/plot_0_Wollastonite.py +++ b/examples/1D_simulation/plot_0_Wollastonite.py @@ -13,6 +13,8 @@ # were obtained from Hansen `et. al.` [#f1]_ import matplotlib as mpl import matplotlib.pyplot as plt +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo from mrsimulator import Simulator from mrsimulator import Site from mrsimulator import SpinSystem @@ -20,8 +22,10 @@ # global plot configuration mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 2 + # %% -# **Step 1** Create the sites. +# **Step 1:** Create the sites. S29_1 = Site( isotope="29Si", isotropic_chemical_shift=-89.0, # in ppm @@ -41,12 +45,12 @@ sites = [S29_1, S29_2, S29_3] # all sites # %% -# **Step 2** Create the spin systems from these sites. Again, we create three +# **Step 2:** Create the spin systems from these sites. Again, we create three # single-site spin systems for better performance. spin_systems = [SpinSystem(sites=[s]) for s in sites] # %% -# **Step 3** Create a Bloch decay spectrum method. +# **Step 3:** Create a Bloch decay spectrum method. method = BlochDecaySpectrum( channels=["29Si"], magnetic_flux_density=14.1, # in T @@ -62,19 +66,32 @@ ) # %% -# **Step 4** Create the Simulator object and add the method and spin-system objects. +# **Step 4:** Create the Simulator object and add the method and spin-system objects. sim_wollastonite = Simulator() sim_wollastonite.spin_systems += spin_systems # add the spin systems sim_wollastonite.methods += [method] # add the method # %% -# **Step 5** Simulate the spectrum. +# **Step 5:** Simulate the spectrum. sim_wollastonite.run() +# The plot of the simulation before post-processing. +ax = plt.subplot(projection="csdm") +ax.plot(sim_wollastonite.methods[0].simulation.real, color="black", linewidth=1) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + # %% -# **Step 6** The plot of the simulation. +# **Step 6:** Add post-simulation processing. +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=70), sp.FFT()] +) +processed_data = post_sim.apply_operations(data=sim_wollastonite.methods[0].simulation) + +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") -ax.plot(sim_wollastonite.methods[0].simulation, color="black", linewidth=1) +ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() diff --git a/examples/1D_simulation/plot_1_PotassiumSulfate.py b/examples/1D_simulation/plot_1_PotassiumSulfate.py index 13e2cf9be..2dd78d419 100644 --- a/examples/1D_simulation/plot_1_PotassiumSulfate.py +++ b/examples/1D_simulation/plot_1_PotassiumSulfate.py @@ -12,6 +12,8 @@ # for :math:`^{33}\text{S}` is obtained from Moudrakovski `et. al.` [#f3]_ import matplotlib as mpl import matplotlib.pyplot as plt +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo from mrsimulator import Simulator from mrsimulator import Site from mrsimulator import SpinSystem @@ -19,8 +21,10 @@ # global plot configuration mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 2 + # %% -# **Step 1** Create the sites, in this case, just the one. +# **Step 1:** Create the sites, in this case, just the one. S33 = Site( name="33S", isotope="33S", @@ -29,11 +33,11 @@ ) # %% -# **Step 2** Create the spin-system from the site. +# **Step 2:** Create the spin-system from the site. spin_system = SpinSystem(sites=[S33]) # %% -# **Step 3** Create a central transition selective Bloch decay spectrum method. +# **Step 3:** Create a central transition selective Bloch decay spectrum method. method = BlochDecayCentralTransitionSpectrum( channels=["33S"], magnetic_flux_density=21.14, # in T @@ -49,19 +53,33 @@ ) # %% -# **Step 4** Create the Simulator object and add the method and the spin-system object. +# **Step 4:** Create the Simulator object and add the method and the spin-system object. sim_K2SO3 = Simulator() sim_K2SO3.spin_systems += [spin_system] # add the spin-system sim_K2SO3.methods += [method] # add the method # %% -# **Step 5** Simulate the spectrum. +# **Step 5:** Simulate the spectrum. sim_K2SO3.run() +# The plot of the simulation before post-processing. +ax = plt.subplot(projection="csdm") +ax.plot(sim_K2SO3.methods[0].simulation.real, color="black", linewidth=1) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + + # %% -# **Step 6** The plot of the simulation. +# **Step 6:** Add post-simulation processing. +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=10), sp.FFT()] +) +processed_data = post_sim.apply_operations(data=sim_K2SO3.methods[0].simulation) + +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") -ax.plot(sim_K2SO3.methods[0].simulation, color="black", linewidth=1) +ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() diff --git a/examples/1D_simulation/plot_2_Coesite.py b/examples/1D_simulation/plot_2_Coesite.py index 68b0fda91..47c4c02f9 100644 --- a/examples/1D_simulation/plot_2_Coesite.py +++ b/examples/1D_simulation/plot_2_Coesite.py @@ -13,6 +13,8 @@ # Grandinetti `et. al.` [#f2]_ import matplotlib as mpl import matplotlib.pyplot as plt +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo from mrsimulator import Simulator from mrsimulator import Site from mrsimulator import SpinSystem @@ -20,8 +22,10 @@ # global plot configuration mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 2 + # %% -# **Step 1** Create the sites. +# **Step 1:** Create the sites. # default unit of isotropic_chemical_shift is ppm and Cq is Hz. O17_1 = Site( @@ -44,14 +48,14 @@ sites = [O17_1, O17_2, O17_3, O17_4, O17_5] # %% -# **Step 2** Create the spin systems from these sites. For optimum performance, we +# **Step 2:** Create the spin systems from these sites. For optimum performance, we # create five single-site spin systems instead of a single five-site spin-system. The # abundance of each spin-system is taken from above reference. abundance = [0.83, 1.05, 2.16, 2.05, 1.90] spin_systems = [SpinSystem(sites=[s], abundance=a) for s, a in zip(sites, abundance)] # %% -# **Step 3** Create a central transition selective Bloch decay spectrum method. +# **Step 3:** Create a central transition selective Bloch decay spectrum method. method = BlochDecayCentralTransitionSpectrum( channels=["17O"], rotor_frequency=14000, # in Hz @@ -71,19 +75,32 @@ # width using 2048 points. # %% -# **Step 4** Create the Simulator object and add the method and spin-system objects. +# **Step 4:** Create the Simulator object and add the method and spin-system objects. sim_coesite = Simulator() sim_coesite.spin_systems += spin_systems # add the spin systems sim_coesite.methods += [method] # add the method # %% -# **Step 5** Simulate the spectrum. +# **Step 5:** Simulate the spectrum. sim_coesite.run() +# The plot of the simulation before post-processing. +ax = plt.subplot(projection="csdm") +ax.plot(sim_coesite.methods[0].simulation.real, color="black", linewidth=1) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + # %% -# **Step 6** The plot of the simulation. +# **Step 6:** Add post-simulation processing. +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=30), apo.Gaussian(sigma=80), sp.FFT()] +) +processed_data = post_sim.apply_operations(data=sim_coesite.methods[0].simulation) + +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") -ax.plot(sim_coesite.methods[0].simulation, color="black", linewidth=1) +ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() diff --git a/examples/Fitting/README.rst b/examples/Fitting/README.rst index 6a1794d59..3b09f909e 100644 --- a/examples/Fitting/README.rst +++ b/examples/Fitting/README.rst @@ -1,3 +1,5 @@ 1D Data Fitting (Least Squares) ------------------------------- -Using LMFIT nonlinear least squares minimization routine to fit a simulation object to experimental data. +The ``mrsimulator`` library is easily integrable with other python-based libraries. +In the following examples, we illustrate the use of LMFIT non-linear least-squares +minimization python package to fit a simulation object to experimental data. diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 77c200f66..dae2855f3 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -1,186 +1,244 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Fitting Cusipidine. -^^^^^^^^^^^^^^^^^^^ +Fitting Cusipidine +================== .. sectionauthor:: Maxwell C. Venetos """ -# sphinx_gallery_thumbnail_number = 3 -#%% -# Often, after obtaining an NMR measurement we must fit tensors to our data so we can -# obtain the tensor parameters. In this example, we will illustrate the use of the *mrsimulator* -# method to simulate the experimental spectrum and fit the simulation to the data allowing us to -# extract the tensor parameters for our spin systems. We will be using the `LMFIT `_ -# methods to establish fitting parameters and fit the spectrum. The following examples will show fitting with -# two synthetic :math:`^{29}\text{Si}` spectra--cuspidine and wollastonite--as well as the -# measurements from an :math:`^{17}\text{O}` experiment on :math:`\text{Na}_{2}\text{SiO}_{3}`. -# The *mrsimulator* library and data make use of CSDM compliant files. -# In this example we will be creating a synthetic spectrum of cuspidine from reported tensor -# parameters and then fit a simulation to the spectrum to demonstrate a simple fitting procedure. -# The :math:`^{29}\text{Si}` tensor parameters were obtained from Hansen `et. al.` [#f1]_ +# %% +# Often, after acquiring an NMR spectrum, we may require some form of least-squares +# analysis to quantify our measurement. A typical recipe for any least-squares analysis +# comprises of two steps: +# +# - Create a "fitting model," and +# - optimize the parameters of the model. +# +# Here, we will use the mrsimulator objects to create a "fitting model," and use the +# `LMFIT `_ library for performing the least-squares +# fitting optimization. +# In this example, we use a synthetic :math:`^{29}\text{Si}` NMR spectrum of cuspidine, +# generated from the tensor parameters reported by Hansen `et. al.` [#f1]_, to +# demonstrate a simple fitting procedure. # # We will begin by importing *matplotlib* and establishing figure size. -import matplotlib.pylab as pylab +import matplotlib as mpl import matplotlib.pyplot as plt -params = {"figure.figsize": (4.5, 3), "font.size": 9} -pylab.rcParams.update(params) +font = {"size": 9} +mpl.rc("font", **font) +mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 3 -#%% -# Next we will import `csdmpy `_ and loading the data file. +# %% +# Import the dataset +# ------------------ +# We will import the `csdmpy `_ +# module and load the synthetic dataset as a CSDM object. import csdmpy as cp filename = "https://osu.box.com/shared/static/a45xj96iekdjrs2beri0nkrow4vjewdh.csdf" synthetic_experiment = cp.load(filename) + +# convert the dimension coordinates from Hz to ppm synthetic_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") -x1, y1 = synthetic_experiment.to_list() +# Normalize the spectrum +synthetic_experiment /= synthetic_experiment.max() -plt.plot(x1, y1) -plt.xlabel("$^{29}$Si frequency / ppm") -plt.xlim(x1.value.max(), x1.value.min()) -plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) +# plot of the synthetic dataset. +ax = plt.subplot(projection="csdm") +ax.plot(synthetic_experiment, color="black", linewidth=1) +ax.set_xlim(-200, 50) +ax.invert_xaxis() plt.tight_layout() plt.show() -#%% -# In order to to fit a simulation to the data we will need to establish a ``Simulator`` object. We will -# use approximate initial parameters to generate our simulation: - -#%% - - -from mrsimulator import Simulator, SpinSystem, Site +# %% +# Create a "fitting model" +# ------------------------ +# +# Before you can fit a simulation to an experiment, in this case, the synthetic dataset, +# you will first need to create a "fitting model." We will use the ``mrsimulator`` +# objects as tools in creating a model for the least-squares fitting. +from mrsimulator import SpinSystem, Simulator from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.post_simulation import PostSimulator -S29 = Site( +# %% +# **Step 1:** Create the guess sites and spin systems. The guess is often based on some +# prior knowledge. For the current example, we know that Cuspidine is a crystalline +# silica polymorph with one crystallographic Si site. Therefore, our guess model is a +# single :math:`^{29}\text{Si}` site spin-system. For non-linear fitting algorithms, +# as a general recommendation, the guess model parameters should be a good starting +# point for the algorithms to converge. + +# the guess model comprising of a single site spin system +site = dict( isotope="29Si", - isotropic_chemical_shift=-80.0, - shielding_symmetric={"zeta": -60, "eta": 0.6}, + isotropic_chemical_shift=-82.0, # in ppm, + shielding_symmetric={"zeta": -63, "eta": 0.4}, # zeta in ppm ) +system_object = SpinSystem( + name="Si Site", + description="A 29Si site in cuspidine", + sites=[site], # from the above code + abundance=100, +) + +# %% +# **Step 2:** Create the method object. The method should be the same as the one used +# in the measurement. In this example, we use the `BlochDecaySpectrum` method. Note, +# when creating the method object, the value of the method parameters must match the +# respective values used in the experiment. method = BlochDecaySpectrum( channels=["29Si"], - magnetic_flux_density=7.1, - rotor_frequency=780, + magnetic_flux_density=7.1, # in T + rotor_frequency=780, # in Hz spectral_dimensions=[ - {"count": 2048, "spectral_width": 25000, "reference_offset": -5000} + { + "count": 2048, + "spectral_width": 25000, # in Hz + "reference_offset": -5000, # in Hz + } ], ) -PS = PostSimulator( - scale=1, apodization=[{"args": [200], "function": "Lorentzian", "dimension": 0}] -) - - +# %% +# **Step 3:** Create the Simulator object and add the method and spin-system objects. sim = Simulator() -sim.spin_systems += [SpinSystem(name="Si29", sites=[S29], abundance=100)] -sim.methods += [method] +sim.spin_systems = [system_object] +sim.methods = [method] sim.methods[0].experiment = synthetic_experiment -sim.methods[0].post_simulation = PS +# %% +# **Step 5:** simulate the spectrum. sim.run() -sim.methods[0].simulation.dimensions[0].to("ppm", "nmr_frequency_ratio") -x, y = sim.methods[0].simulation.to_list() -y = sim.methods[0].apodize().real +# %% +# **Step 6:** Create a SignalProcessor class and apply post simulation processing. +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo + +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=200), sp.FFT(), sp.Scale(factor=1.5)] +) +processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) -plt.plot(x, y) -plt.xlabel("$^{29}$Si frequency / ppm") -plt.xlim(x.value.max(), x.value.min()) -plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) +# %% +# **Step 7:** The plot the spectrum. We also plot the synthetic dataset for comparison. +ax = plt.subplot(projection="csdm") +ax.plot(processed_data.real, c="k", linewidth=1, label="guess spectrum") +ax.plot(synthetic_experiment.real, c="r", linewidth=1.5, alpha=0.5, label="experiment") +ax.set_xlim(-200, 50) +ax.invert_xaxis() +plt.legend() plt.tight_layout() plt.show() -#%% -# Next, we will need a list of parameters that will be used in the fit. the *LMFIT* library allows us to create -# a list of parameters rather easily using the `Parameters `_ class. -# We have created a function to parse the -# ``simulator`` object for available parameters and construct an *LMFIT* ``Parameter`` object which is shown in the next two -# examples on fitting. Here, however, we will construct the parameter list explicitly to demonstrate how the parameters -# are created. - -#%% - - -from lmfit import Minimizer, Parameters +# %% +# Setup a Least-squares minimization +# ---------------------------------- +# +# Now that our model is ready, the next step is to set up a least-squares minimization. +# You may use any optimization package of choice, here we show an application using +# LMFIT. You may read more on the LMFIT +# `documentation page `_. +# +# Create fitting parameters +# ''''''''''''''''''''''''' +# +# Next, you will need a list of parameters that will be used in the fit. The *LMFIT* +# library provides a `Parameters `_ +# class to create a list of parameters rather easily. +from lmfit import Minimizer, Parameters, fit_report +site1 = system_object.sites[0] params = Parameters() -params.add(name="iso", value=-80) -params.add(name="eta", value=0.6, min=0, max=1) -params.add(name="zeta", value=-60) -params.add(name="sigma", value=200) -params.add(name="scale", value=1) -params +params.add(name="iso", value=site1.isotropic_chemical_shift) +params.add(name="eta", value=site1.shielding_symmetric.eta, min=0, max=1) +params.add(name="zeta", value=site1.shielding_symmetric.zeta) +params.add(name="FWHM", value=post_sim.operations[1].FWHM) +params.add(name="factor", value=post_sim.operations[3].factor) -#%% -# We will next set up an error function that will update the simulation throughout the minimization. -# We will construct a simple function here to demonstrate the *LMFIT* library, however, the next examples -# will showcase a fitting function provided in the *mrsimulator* library which automates the process. +# %% +# Create a minimization function +# '''''''''''''''''''''''''''''' +# +# Note, the above set of parameters are does not know the model. You will need to set up +# a function that will +# +# - update the parameters of the `Simulator` and `SignalProcessor` object based on the +# LMFIT parameter updates, +# - re-simulate the spectrum based on the updated values, and +# - return the difference between the experiment and simulation. -def test_function(params, data, sim): +def minimization_function(params, sim, post_sim): values = params.valuesdict() - intensity = data.dependent_variables[0].components[0].real + # the experiment data as a Numpy array + intensity = sim.methods[0].experiment.dependent_variables[0].components[0].real # Here, we update simulation parameters iso, eta, and zeta for the site object site = sim.spin_systems[0].sites[0] site.isotropic_chemical_shift = values["iso"] site.shielding_symmetric.eta = values["eta"] site.shielding_symmetric.zeta = values["zeta"] - sim.methods[0].post_simulation.scale = values["scale"] - sim.methods[0].post_simulation.apodization[0].args = [values["sigma"]] - # here we run the simulation + # run the simulation sim.run() - # here we apodize the signal to simulate line broadening - y = sim.methods[0].apodize().real - - y_factored = y * values["scale"] - - return intensity - y_factored + # update the SignalProcessor parameter and apply line broadening. + post_sim.operations[3].factor = values["factor"] + post_sim.operations[1].FWHM = values["FWHM"] + processed_data = post_sim.apply_operations(sim.methods[0].simulation) -#%% -# With the synthetic data, simulation, and the parameters we are ready to perform the fit. To fit, we use -# the *LMFIT* `Minimizer `_ class. + # return the difference vector. + return intensity - processed_data.dependent_variables[0].components[0].real -#%% -minner = Minimizer(test_function, params, fcn_args=(synthetic_experiment, sim)) +# %% +# .. note:: +# To automate the fitting process, we provide a function to parse +# the ``simulator`` object for parameters and construct an *LMFIT* ``Parameters`` +# object. Similarly, a minimization function, analogous to the above +# `minimization_function`, is also included in the *mrsimulator* library. See the +# next example for usage instructions. +# +# Perform the least-squares minimization +# '''''''''''''''''''''''''''''''''''''' +# +# With the synthetic data, simulation, and the parameters, we are ready to perform the +# fit. To fit, we use the *LMFIT* +# `Minimizer `_ class. One consideration +# for the case of the magic-angle spinning fitting is we must use a discrete +# minimization method, such as 'powell', as the chemical shift varies discretely. +minner = Minimizer(minimization_function, params, fcn_args=(sim, post_sim)) result = minner.minimize(method="powell") -result - -#%% -# After the fit, we can plot the new simulated spectrum. - +print(fit_report(result)) +# %% +# The plot of the fit, measurement and the residuals is shown below. plt.figsize = (4, 3) -residual = synthetic_experiment.copy() -residual[:] = result.residual -plt.plot(*synthetic_experiment.to_list(), label="Spectrum") -plt.plot(*(synthetic_experiment - residual).to_list(), "r", alpha=0.5, label="Fit") -plt.plot(*residual.to_list(), alpha=0.5, label="Residual") +x, y_data = synthetic_experiment.to_list() +residual = result.residual +plt.plot(x, y_data, label="Spectrum") +plt.plot(x, y_data - residual, "r", alpha=0.5, label="Fit") +plt.plot(x, residual, alpha=0.5, label="Residual") -plt.xlim(x1.value.max(), x1.value.min()) plt.xlabel("Frequency / Hz") +plt.xlim(-200, 50) +plt.gca().invert_xaxis() plt.grid(which="major", axis="both", linestyle="--") plt.legend() - plt.tight_layout() plt.show() -#%% -# +# %% # .. [#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/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 1ddb2e563..f0d4b7d70 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -1,178 +1,214 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Fitting Sodium Silicate. -^^^^^^^^^^^^^^^^^^^^^^^^ +Fitting Crystalline Sodium Metasilicate +======================================= .. sectionauthor:: Maxwell C. Venetos """ -# sphinx_gallery_thumbnail_number = 3 -#%% -# Often, after obtaining an NMR measurement we must fit tensors to our data so we can -# obtain the tensor parameters. In this example, we will illustrate the use of the *mrsimulator* -# method to simulate the experimental spectrum and fit the simulation to the data allowing us to -# extract the tensor parameters for our spin systems. We will be using the `LMFIT `_ -# methods to establish fitting parameters and fit the spectrum. The following examples will show fitting with -# two synthetic :math:`^{29}\text{Si}` spectra--cuspidine and wollastonite--as well as the -# measurements from an :math:`^{17}\text{O}` experiment on :math:`\text{Na}_{2}\text{SiO}_{3}`. -# The *mrsimulator* library and data make use of CSDM compliant files. -# In this example we will fit a simulation to an experimentally obtained :math:`^{17}\text{O}` spectrum. -# We use the :math:`^{17}\text{O}` tensor information from Grandinetti `et. al.` [#f5]_ +# %% +# In this example, we illustrate the use of the mrsimulator objects to +# +# - create a spin-system model, +# - use the model to perform a least-squares fit on the experimental, and +# - extract the tensor parameters of the spin - system model. +# +# We will be using the `LMFIT `_ methods to +# establish fitting parameters and fit the spectrum. The following example illustrates +# the least-squares fitting on a :math:`^{17}\text{O}` measurement of +# :math:`\text{Na}_{2}\text{SiO}_{3}` [#f5]_. # -# We will begin by importing *matplotlib* and establishing figure size. -import matplotlib.pylab as pylab +# We will begin by importing relevant modules and presetting figure style and layout. +import csdmpy as cp +import matplotlib as mpl import matplotlib.pyplot as plt +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo +from lmfit import Minimizer +from lmfit import report_fit +from mrsimulator import Simulator +from mrsimulator import Site +from mrsimulator import SpinSystem +from mrsimulator.methods import BlochDecayCentralTransitionSpectrum +from mrsimulator.spectral_fitting import LMFIT_min_function +from mrsimulator.spectral_fitting import make_LMFIT_parameters -params = {"figure.figsize": (4.5, 3), "font.size": 9} -pylab.rcParams.update(params) +font = {"size": 9} +mpl.rc("font", **font) +mpl.rcParams["figure.figsize"] = [4.25, 3.0] +# sphinx_gallery_thumbnail_number = 3 -#%% -# Next we will import `csdmpy `_ and loading the data file. -import csdmpy as cp +# %% +# Import the dataset +# ------------------ +# +# Import the experimental data. In this example, we will import the data file serialized +# with the CSDM file-format, using the +# `csdmpy `_ module. +filename = "https://osu.box.com/shared/static/kfgt0jxgy93srsye9pofdnoha6qy58qf.csdf" +oxygen_experiment = cp.load(filename) +# For spectral fitting, we only focus on the real part of the complex dataset +oxygen_experiment = oxygen_experiment.real -filename = "https://osu.box.com/shared/static/kfgt0jxgy93srsye9pofdnoha6qy58qf.csdf" -oxygen_experiment = cp.load(filename).real +# Convert the dimension coordinates from Hz to ppm. oxygen_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") -x1, y1 = oxygen_experiment.to_list() +# Normalize the spectrum +oxygen_experiment /= oxygen_experiment.max() -plt.plot(x1, y1) -plt.xlabel("$^{17}$O frequency / ppm") -plt.xlim(x1.value.max(), x1.value.min()) -plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) +# plot of the dataset. +ax = plt.subplot(projection="csdm") +ax.plot(oxygen_experiment, color="black", linewidth=1) +ax.set_xlim(-50, 100) +ax.invert_xaxis() plt.tight_layout() plt.show() - -#%% -# Next, we will want to create a ``simulator`` object that we will use to fit to our -# spectrum. We will need to import the necessary libraries for the *mrsimulator* -# methods. We will then create ``SpinSystem`` objects. - -#%% - -from mrsimulator import Simulator, SpinSystem, Site -from mrsimulator.methods import BlochDecayCentralTransitionSpectrum -from mrsimulator.post_simulation import PostSimulator - -sim = Simulator() +# %% +# Create a "fitting model" +# ------------------------ +# +# Next, we will create a ``simulator`` object that we use to fit the spectrum. We will +# start by creating the guess ``SpinSystem`` objects. +# +# **Step 1:** Create the guess sites. O17_1 = Site( isotope="17O", - isotropic_chemical_shift=60.0, - quadrupolar={"Cq": 4200000, "eta": 0.5}, + isotropic_chemical_shift=60.0, # in ppm, + quadrupolar={"Cq": 4.2e6, "eta": 0.5}, # Cq in Hz ) + O17_2 = Site( - isotope="17O", isotropic_chemical_shift=40, quadrupolar={"Cq": 2400000, "eta": 0.18} + isotope="17O", + isotropic_chemical_shift=40.0, # in ppm, + quadrupolar={"Cq": 2.4e6, "eta": 0.5}, # Cq in Hz ) -spin_systems = [SpinSystem(sites=[site]) for site in [O17_1, O17_2]] - - +# %% +# **Step 2:** Create the spin-systems for the guess sites. +system_object = [SpinSystem(sites=[s]) for s in [O17_1, O17_2]] # from the above code + +# %% +# **Step 3:** Create the method object. Note, when performing the least-squares fit, you +# must create an appropriate method object which matches the method used in acquiring +# the experimental data. The attribute values of this method must be set to match the +# exact conditions under which the experiment was acquired. This including the +# acquisition channels, the magnetic flux density, rotor angle, rotor frequency, and +# the spectral/spectroscopic dimension. In the following example, we set up a central +# transition selective Bloch decay spectrum method, where the spectral/spectroscopic +# information is obtained from the metadata of the CSDM dimension. The remaining +# attribute values are set to the experimental conditions. + +# get the count, increment, and coordinates_offset info from the experiment dimension. count = oxygen_experiment.dimensions[0].count increment = oxygen_experiment.dimensions[0].increment.to("Hz").value offset = oxygen_experiment.dimensions[0].coordinates_offset.to("Hz").value method = BlochDecayCentralTransitionSpectrum( channels=["17O"], - magnetic_flux_density=9.4, - rotor_frequency=14000, + magnetic_flux_density=9.4, # in T + rotor_frequency=14000, # in Hz spectral_dimensions=[ { "count": count, - "spectral_width": count * increment, - "reference_offset": offset, + "spectral_width": count * increment, # in Hz + "reference_offset": offset, # in Hz } ], ) +# %% +# Now we add the experimental data to the ``experiment`` attribute of the above method. +method.experiment = oxygen_experiment -PS = PostSimulator( - scale=oxygen_experiment.dependent_variables[0].components[0].max().real / 4, - apodization=[{"args": [100], "function": "Lorentzian", "dimension": 0}], -) - -sim.spin_systems += spin_systems -sim.methods += [method] -sim.methods[0].post_simulation = PS -sim.methods[0].experiment = oxygen_experiment +# %% +# **Step 4:** Create the Simulator object and add the method and spin-system objects. +sim = Simulator() +sim.spin_systems = system_object +sim.methods = [method] -# To avoid querying at every iteration we will save the relevant transition pathways -for iso in spin_systems: +# %% +# **Step 5:** Simulate the spectrum. +for iso in sim.spin_systems: + # To avoid querying at every iteration, save the relevant transition pathways info iso.transition_pathways = method.get_transition_pathways(iso).tolist() -sim.run() - -sim.methods[0].simulation.dimensions[0].to("ppm", "nmr_frequency_ratio") -x, y = sim.methods[0].simulation.to_list() -y = sim.methods[0].apodize().real +sim.run() -plt.plot(x, y) -plt.xlabel("$^{17}$O frequency / ppm") -plt.xlim(x.value.max(), x.value.min()) -plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) +# %% +# **Step 6:** Create the SignalProcessor class object and add and apply the list of +# post-simulation operations. +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=40), sp.FFT(), sp.Scale(factor=0.6)] +) +processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) + +# %% +# **Step 7:** The plot the guess simulation (black) along with the spectrum (red) is +# shown below. +ax = plt.subplot(projection="csdm") +ax.plot(processed_data.real, color="black", linewidth=1, label="guess spectrum") +ax.plot(oxygen_experiment.real, c="r", linewidth=1.5, alpha=0.5, label="experiment") +ax.set_xlim(-50, 100) +ax.invert_xaxis() +plt.legend() plt.tight_layout() plt.show() -#%% -# Once we have our simulation we must create our list of parameters to use in our -# fitting. We will be using the `Parameters `_ class from *LMFIT*. +# %% +# Least-squares minimization with LMFIT +# ------------------------------------- # -# In each experiment the number of spin systems and sites present as well as their -# attributes may vary. To simplify the parameter list creation we will use the -# :func:`~mrsimulator.spectral_fitting.make_fitting_parameters` - -#%% - - -from mrsimulator.spectral_fitting import make_fitting_parameters - -params = make_fitting_parameters(sim) -params - -#%% -# With an experimental spectrum, a simulaton, and a list of parameters we are now -# ready to perform a fit. This fit will be performed using the *LMFIT* library as -# well as our error function, :func:`~mrsimulator.spectral_fitting.min_function`. The arguments -# for ``min_function`` are the intensities from the experimental data and the simulation -# CSDM object. Reporting the results of the fit will give us our tensor parameters. +# Once you have a fitting model, you need to create the list of parameters to use in the +# least-squares fitting. For this, you may use the +# `Parameters `_ class from *LMFIT*, +# as described in the previous example. +# Here, we make use of a utility function, +# :func:`~mrsimulator.spectral_fitting.make_LMFIT_parameters`, that considerably +# simplifies the generation of the parameters object. # -# One thing to note is that the names of our parameters must correspond to their addresses within the simulation object -# in order to update the simulation during the fit. The *LMFIT* library does not allow for the use of special characters -# such as "\[", "\]", or "." so our current workaround is converting the special characters to their corresponding HTML -# character code numbers and converting back to the special character when updating the simulation. - - -#%% - -from mrsimulator.spectral_fitting import min_function -from lmfit import Minimizer +# **Step 8:** Create a list of parameters. +params = make_LMFIT_parameters(sim, post_sim) -minner = Minimizer(min_function, params, fcn_args=(sim, "Lorentzian")) +# %% +# The `make_LMFIT_parameters` parses the instances of the ``Simulator`` and the +# ``PostSimulator`` objects for parameters and returns an LMFIT `Parameters` object. +# +# **Customize the Parameters:** +# You may customize the parameters list, ``params``, as desired. Here, we add a +# constraint on the fit by fixing the site abundances for the spin-systems at index +# 1 and 2 to 50%. +params["sys_0_abundance"].vary = False # fix the abundance +params.pretty_print() + +# %% +# **Step 9:** Perform least-squares minimization. For the user's convenience, we also +# provide a utility function, :func:`~mrsimulator.spectral_fitting.LMFIT_min_function`, +# for evaluating the difference vector between the simulation and experiment, based on +# the parameters update. You may use this function directly as the argument of the +# LMFIT Minimizer class, as follows, +minner = Minimizer(LMFIT_min_function, params, fcn_args=(sim, post_sim)) result = minner.minimize() -result - -#%% -# Next, we can compare the fit to the experimental data: - -#%% - +report_fit(result) +# %% +# **Step 10:** The plot of the fit, measurement and the residuals is shown below. plt.figsize = (4, 3) -residual = oxygen_experiment.copy() -residual[:] = result.residual -plt.plot(*oxygen_experiment.to_list(), label="Spectrum") -plt.plot(*(oxygen_experiment - residual).to_list(), "r", alpha=0.5, label="Fit") -plt.plot(*residual.to_list(), alpha=0.5, label="Residual") +x, y_data = oxygen_experiment.to_list() +residual = result.residual +plt.plot(x, y_data, label="Spectrum") +plt.plot(x, y_data - residual, "r", alpha=0.5, label="Fit") +plt.plot(x, residual, alpha=0.5, label="Residual") -plt.xlim(x.value.max(), x.value.min()) plt.xlabel("$^{17}$O frequency / ppm") +plt.xlim(-50, 100) +plt.gca().invert_xaxis() plt.grid(which="major", axis="both", linestyle="--") plt.legend() - plt.tight_layout() plt.show() -#%% +# %% # # .. [#f5] T. M. Clark, P. Florian, J. F. Stebbins, and P. J. Grandinetti, # An :math:`^{17}\text{O}` NMR Investigation of Crystalline Sodium Metasilicate: diff --git a/requirements-dev.txt b/requirements-dev.txt index 32095f9a5..b94955b8a 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,6 +1,6 @@ numpy>=1.17 matplotlib>=3.0 -csdmpy>=0.3.0 +csdmpy>=0.3.1 pydantic>=1.0 monty>=2.0.4 typing-extensions>=3.7 diff --git a/requirements.txt b/requirements.txt index 15cfc400b..660ae1672 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy>=1.17 matplotlib>=3.0 -csdmpy>=0.3.0 +csdmpy>=0.3.1 pydantic>=1.0 monty>=2.0.4 typing-extensions>=3.7 diff --git a/setup.py b/setup.py index cb99984f3..f3cd78984 100644 --- a/setup.py +++ b/setup.py @@ -257,7 +257,7 @@ def message(lib): setup_requires=["numpy>=1.17"], install_requires=[ "numpy>=1.17", - "csdmpy>=0.3", + "csdmpy>=0.3.1", "pydantic>=1.0", "monty>=2.0.4", "typing-extensions>=3.7", diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 4c183fe1f..b30933370 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -7,7 +7,6 @@ import csdmpy as cp import numpy as np -from mrsimulator.post_simulation import PostSimulator from mrsimulator.spin_system.isotope import Isotope from mrsimulator.transition import Transition from mrsimulator.transition.transition_list import TransitionList @@ -103,14 +102,12 @@ class Method(Parseable): >>> bloch.description 'Huh!' - post_simulation: An optional dict with post-simulation parameters. """ name: str = None label: str = None description: str = None channels: List[str] = [] spectral_dimensions: List[SpectralDimension] = [{}] - post_simulation: PostSimulator = None simulation: Union[cp.CSDM, np.ndarray] = None experiment: Union[cp.CSDM, np.ndarray] = None @@ -208,8 +205,6 @@ def to_dict_with_units(self): def dict(self, **kwargs): temp_dict = super().dict(**kwargs) - if self.post_simulation is not None: - temp_dict["post_simulation"] = self.post_simulation.dict() if self.simulation is not None: temp_dict["simulation"] = self.simulation.to_dict(update_timestamp=True) if self.experiment is not None: @@ -331,28 +326,3 @@ def get_transition_pathways(self, spin_system) -> np.ndarray: # segments.append(np.asarray(P_segment)) # append the intersection # return cartesian_product(*segments) - - def apodize(self, **kwargs): - """Returns the result of passing the selected apodization function . - - Args: - csdm: simulation object - - Returns: - A Numpy array. - """ - - csdm = self.simulation - for dim in csdm.dimensions: - dim.to("Hz", "nmr_frequency_ratio") - apo = self.post_simulation.apodization - - sum_ = 0 - - for item in apo: - sum_ += item.apodize(csdm) - - for dim in csdm.dimensions: - dim.to("ppm", "nmr_frequency_ratio") - - return self.post_simulation.scale * sum_ diff --git a/src/mrsimulator/post_simulation.py b/src/mrsimulator/post_simulation.py deleted file mode 100644 index a8f8d1e33..000000000 --- a/src/mrsimulator/post_simulation.py +++ /dev/null @@ -1,98 +0,0 @@ -# -*- coding: utf-8 -*- -"""The Event class.""" -from typing import List -from typing import Optional - -import numpy as np -from mrsimulator.utils.parseable import Parseable -from numpy.fft import fft -from numpy.fft import fftshift -from numpy.fft import ifft -from numpy.fft import ifftshift - - -class Apodization(Parseable): - """The Apodization class.""" - - function: str - dimension: int = 0 - args: List[float] = [0] - fraction: Optional[float] = 1 - - class Config: - validate_assignment = True - - def apodize(self, csdm, **kwargs): - """Returns the result of passing the selected apodization function . - - Args: - csdm: simulation object - - Returns: - A Numpy array - """ - function_mapping = {"Gaussian": Gaussian, "Lorentzian": Lorentzian} - - y = csdm.dependent_variables[0].components[0] - # x = csdm.dimensions[self.dimension].coordinates - axis = -self.dimension - 1 - fapp = function_mapping[self.function]( - csdm, self.dimension, self.args, **kwargs - ) - - recip_dim = csdm.dimensions[0].reciprocal - coord = csdm.dimensions[0].coordinates - phase = np.exp(2j * np.pi * recip_dim.coordinates_offset * coord) - - TimeDomain = ifft(ifftshift(y * phase, axes=axis), axis=axis) - - appodized = TimeDomain * fapp - fft_output = fftshift(fft(appodized, axis=axis), axes=axis).real - - # csdm.dependent_variables[index].components[0] = ( - # self.fraction * phase.conj() * fft_output - # ) - - return self.fraction * phase.conj() * fft_output - - -def Lorentzian(csdm, axis, arg, **kwargs): - """Lorentzian apodization function. - - Args: - Self: simulation object - arg: list with one entry. The full-width-half-max in Hz of the Lorentzian function - - Returns: - A Numpy array - """ - inv_x = csdm.dimensions[axis].reciprocal_coordinates().to("s").value - return np.exp(-arg[0] * np.pi * np.abs(inv_x)) - - -def Gaussian(csdm, axis, arg, **kwargs): - """Gaussian apodization function. - - Args: - Self: simulation object - sigma: list with one entry. standard deviation in Hz of the Gaussian function - - Returns: - A Numpy array - """ - inv_x = csdm.dimensions[axis].reciprocal_coordinates().to("s").value - return np.exp(-((inv_x * arg[0] * np.pi) ** 2) * 2) - - -class PostSimulator(Parseable): - """scale: scaling factor - """ - - scale: Optional[float] = 1.0 - apodization: List[Apodization] = [] - truncation_artifact: Optional[bool] = False - dead_time: Optional[float] = 0.0 - - class Config: - validate_assignment = True - arbitrary_types_allowed = True diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py new file mode 100644 index 000000000..df3250531 --- /dev/null +++ b/src/mrsimulator/signal_processing/__init__.py @@ -0,0 +1,168 @@ +# -*- coding: utf-8 -*- +"""The Event class.""" +from sys import modules +from typing import List + +import csdmpy as cp +from pydantic import BaseModel + +from ._base import AbstractOperation + +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" + + +class SignalProcessor(BaseModel): + """ + Signal processing class to apply a series of operations to the dependent variables + of the simulation dataset. + + Attributes + ---------- + + operations: List + A list of operations. + + Examples + -------- + + >>> post_sim = SignalProcessor(operations=[o1, o2]) # doctest: +SKIP + """ + + processed_data: cp.CSDM = None + operations: List[AbstractOperation] = [] + + class Config: + validate_assignment = True + arbitrary_types_allowed = True + + @classmethod + def parse_dict_with_units(self, py_dict): + """Parse a list of operations dictionary to a SignalProcessor class object. + + Args: + pt_dict: A python dict object. + """ + lst = [] + for op in py_dict["operations"]: + if op["function"] == "apodization": + lst.append( + getattr( + modules[__name__].apodization, op["type"] + ).parse_dict_with_units(op) + ) + else: + lst.append( + getattr(modules[__name__], op["function"]).parse_dict_with_units(op) + ) + return SignalProcessor(operations=lst) + + def to_dict_with_units(self): + """ + Serialize the SignalProcessor object to a JSON compliant python dictionary + object, where physical quantities are represented as string with a value and a + unit. + + Returns: + A Dict object. + """ + lst = [] + for i in self.operations: + lst += [i.to_dict_with_units()] + op = {} + + op["operations"] = lst + return op + + def apply_operations(self, data, **kwargs): + """ + Function to apply all the operation functions in the operations member of a + SignalProcessor object. Operations applied sequentially over the data member. + + Returns: + CSDM object: A copy of the data member with the operations applied to it. + """ + if not isinstance(data, cp.CSDM): + raise ValueError("The data must be a CSDM object.") + copy_data = data.copy() + for filters in self.operations: + copy_data = filters.operate(copy_data) + self.processed_data = copy_data + + return copy_data + + +class Scale(AbstractOperation): + """ + Scale the amplitudes of the dependent variables from the simulation data object. + + Args: + float factor: The scaling factor applied to all the dependent variables in the + simulation data. The default value is 1. + + Example + ------- + + >>> import mrsimulator.signal_processing as sp + >>> operation1 = sp.Scale(factor=20) + """ + + factor: float = 1 + + def operate(self, data): + r"""Applies the operation for which the class is named for. + + .. math:: + f(\vec(x)) = scale*\vec(x) + + Args: + data: CSDM object + """ + data *= self.factor + return data + + +class IFFT(AbstractOperation): + """ + Apply an inverse Fourier transform on the dependent variables from the simulation + data object. + + Args: + int dim_indx: Data dimension index along which the function is applied. + + Example + ------- + + >>> operation2 = sp.IFFT(dim_indx=0) + """ + + dim_indx: int = 0 + + def operate(self, data): + """Applies the operation for which the class is named for. + + Args: + data: CSDM object + """ + return data.fft(axis=self.dim_indx) + + +class FFT(IFFT): + """ + Apply a forward Fourier transform on the dependent variables from the simulation + data object. + + Args: + int dim_indx: Data dimension index along which the function is applied. + + Example + ------- + + >>> operation3 = sp.FFT(dim_indx=0) + """ + + pass + + +class complex_conjugate(AbstractOperation): + pass diff --git a/src/mrsimulator/signal_processing/_base.py b/src/mrsimulator/signal_processing/_base.py new file mode 100644 index 000000000..d3dd21287 --- /dev/null +++ b/src/mrsimulator/signal_processing/_base.py @@ -0,0 +1,34 @@ +# -*- coding: utf-8 -*- +"""The AbstractOperation class.""" +from mrsimulator.utils.parseable import Parseable + + +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" + + +class AbstractOperation(Parseable): + """A base class for signal processing operations.""" + + @property + def function(self): + return self.__class__.__name__ + + def to_dict_with_units(self): + my_dict = super().to_dict_with_units() + my_dict["function"] = self.function + if hasattr(self, "type"): + my_dict["type"] = self.type + return my_dict + + @classmethod + def parse_dict_with_units(cls, py_dict): + """Parse dictionary for SignalProcessor + + Args: + py_dict (dict): A python dictionary representation of the operation. + """ + my_dict_copy = py_dict.copy() + if "function" in my_dict_copy.keys(): + my_dict_copy.pop("function") + return super().parse_dict_with_units(my_dict_copy) diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py new file mode 100644 index 000000000..20392b632 --- /dev/null +++ b/src/mrsimulator/signal_processing/apodization.py @@ -0,0 +1,211 @@ +# -*- coding: utf-8 -*- +"""The Event class.""" +from sys import modules +from typing import ClassVar +from typing import Dict +from typing import Union + +import numpy as np + +from ._base import AbstractOperation + +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" + + +class AbstractApodization(AbstractOperation): + dim_indx: int = 0 + dep_var_indx: Union[int, list, tuple] = None # if none apply to all + + @classmethod + def parse_dict_with_units(cls, py_dict): + obj = super().parse_dict_with_units(py_dict) + return getattr(modules[__name__], py_dict["type"])(**obj.dict()) + + @property + def function(self): + return "apodization" + + @property + def type(self): + """The type apodization function.""" + return self.__class__.__name__ + + def set_property_units(self, unit, prop): + """ + Populate the property unit attribute of the class based on the dimension unit. + """ + if unit.physical_type == "frequency": + self.property_units = {f"{prop}": "s"} + return + self.property_units = {f"{prop}": "Hz"} + + @staticmethod + def _get_dv_indexes(indexes, n): + """Return a list of dependent variable indexes. + + Args: + indexes: An interger, list of integers, or None indicating the dv indexes. + n: Total number of dependent variables in the CSDM object. + """ + if indexes is None: + return np.arange(n) + if isinstance(indexes, int): + return [indexes] + if isinstance(indexes, (list, tuple)): + return np.asarray(indexes) + + def _operate(self, data, fn, prop_name, prop_value): + """A generic operation function. + + Args: + data: A CSDM object. + fn: The apodization function. + prop_name: The argument name for the function fn. + prop_value: The argument value for the function fn. + """ + x = data.dimensions[self.dim_indx].coordinates + x_value, x_unit = x.value, x.unit + + self.set_property_units(x_unit, prop_name) + + apodization_vactor = fn(x_value, prop_value) + + n = len(data.dependent_variables) + dv_indexes = self._get_dv_indexes(self.dep_var_indx, n=n) + + for i in dv_indexes: + data.dependent_variables[i].components[0] *= apodization_vactor + return data + + +class Gaussian(AbstractApodization): + r"""Apodize a dependent variable of the simulation data object by a Gaussian + function. The function follows + + .. math:: + f(x) = e^{-2 x^2 \sigma^2 \pi^2}, + + where :math:`x` are the coordinates of the data dimension and :math:`\sigma` is + the standard deviation. + + Args: + int dim_indx: Data dimension index to apply the function along. The default + value is 0. + float sigma: The standard deviation, :math:`\sigma`, of the Gaussian function. + The default value is 0. + int dep_var_indx: Data dependent variable index to apply the function to. If + the type None, the operation will be applied to every dependent variable. + + Example + ------- + + >>> import mrsimulator.signal_processing.apodization as apo + >>> operation4 = apo.Gaussian(sigma=143.4, dim_indx=0, dep_var_indx=0) + """ + + sigma: float = 0 + + property_unit_types: ClassVar = {"sigma": ["time", "frequency"]} + property_default_units: ClassVar = {"sigma": ["s", "Hz"]} + property_units: Dict = {"sigma": "Hz"} + + @staticmethod + def fn(x, arg): + return np.exp(-2 * ((x * arg * np.pi) ** 2)) + + def operate(self, data): + """ + Applies the operation for which the class is named for. + + data: CSDM object + dep_var: int. The index of the dependent variable to apply operation to + """ + + return self._operate(data, fn=self.fn, prop_name="sigma", prop_value=self.sigma) + + +# class ExponentialAbs(AbstractApodization): +# r"""Apodize a dependent variable of the simulation data object by an exponential +# function. The function follows + +# .. math:: +# f(x) = e^{-\Tau |x| \pi}, + +# where :math:`x` are the coordinates of the data dimension and :math:`\Tau` is +# the width parameter. + +# Args: +# int dim_indx: Data dimension index to apply the function along. +# float FWHM: The full width at half maximum parameter, :math:`\Tau`. +# int dep_var_indx: Data dependent variable index to apply the function to. If +# the type None, the operation will be applied to every dependent variable. + +# Example +# ------- + +# >>> operation5 = apo.Exponential(sigma=143.4, dim_indx=0, dep_var_indx=0) +# """ + +# FWHM: float = 0 + +# property_unit_types: ClassVar = {"FWHM": ["time", "frequency"]} +# property_default_units: ClassVar = {"FWHM": ["s", "Hz"]} +# property_units: Dict = {"FWHM": "Hz"} + +# @staticmethod +# def fn(x, arg): +# return np.exp(-arg * np.pi * np.abs(x)) + +# def operate(self, data): +# """ +# Applies the operation for which the class is named for. + +# data: CSDM object +# dep_var: int. The index of the dependent variable to apply operation to +# """ + +# return self._operate(data, fn=self.fn, prop_name="FWHM", prop_value=self.FWHM) + + +class Exponential(AbstractApodization): + r"""Apodize a dependent variable of the simulation data object by an exponential + function. The function follows + + .. math:: + f(x) = e^{-\Tau |x| \pi}, + + where :math:`x` are the coordinates of the data dimension and :math:`\Tau` is + the width parameter. + + Args: + int dim_indx: Data dimension index to apply the function along. + float FWHM: The full width at half maximum parameter, :math:`\Tau`. + int dep_var_indx: Data dependent variable index to apply the function to. If + the type None, the operation will be applied to every dependent variable. + + Example + ------- + + >>> operation5 = apo.Exponential(FWHM=143.4, dim_indx=0, dep_var_indx=0) + """ + + FWHM: float = 0 + + property_unit_types: ClassVar = {"FWHM": ["time", "frequency"]} + property_default_units: ClassVar = {"FWHM": ["s", "Hz"]} + property_units: Dict = {"FWHM": "Hz"} + + @staticmethod + def fn(x, arg): + return np.exp(-arg * np.pi * np.abs(x)) + + def operate(self, data): + """ + Applies the operation for which the class is named for. + + data: CSDM object + dep_var: int. The index of the dependent variable to apply operation to + """ + + return self._operate(data, fn=self.fn, prop_name="FWHM", prop_value=self.FWHM) diff --git a/src/mrsimulator/signal_processing/tests/test_apodization.py b/src/mrsimulator/signal_processing/tests/test_apodization.py new file mode 100644 index 000000000..db6792af3 --- /dev/null +++ b/src/mrsimulator/signal_processing/tests/test_apodization.py @@ -0,0 +1,178 @@ +# -*- coding: utf-8 -*- +"""Apodization test""" +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo +import numpy as np +from mrsimulator import Simulator +from mrsimulator import SpinSystem +from mrsimulator.methods import BlochDecaySpectrum + +sim = Simulator() +the_site = {"isotope": "1H", "isotropic_chemical_shift": "0 ppm"} +the_spin_system = {"name": "site A", "sites": [the_site], "abundance": "80%"} +spin_system_object = SpinSystem.parse_dict_with_units(the_spin_system) +sim.spin_systems += [spin_system_object, spin_system_object, spin_system_object] +sim.config.decompose_spectrum = "spin_system" + +method_1 = BlochDecaySpectrum( + channels=["1H"], + magnetic_flux_density=9.4, + rotor_angle=0, + rotor_frequency=0, + spectral_dimensions=[ + {"count": 4096, "spectral_width": 25000, "reference_offset": 0} + ], +) + +PS_0 = [sp.Scale(factor=10)] + +PS_1 = [ + sp.IFFT(dim_indx=0), + apo.Exponential(FWHM=200, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), +] + +PS_2 = [ + sp.IFFT(dim_indx=0), + apo.Gaussian(sigma=20, dim_indx=0, dep_var_indx=[0, 1]), + sp.FFT(dim_indx=0), +] + +PS_3 = [ + sp.IFFT(dim_indx=0), + apo.Gaussian(sigma=20, dim_indx=0, dep_var_indx=None), + sp.FFT(dim_indx=0), +] + +sim.methods += [method_1] +sim.run() + + +freqHz = sim.methods[0].spectral_dimensions[0].coordinates_Hz() + + +def test_scale(): + post_sim = sp.SignalProcessor(operations=PS_0) + data = post_sim.apply_operations(data=sim.methods[0].simulation) + _, y0, y1, y2 = sim.methods[0].simulation.to_list() + _, y0_, y1_, y2_ = data.to_list() + + assert y0_.max() / y0.max() == 10, "Scaling failed" + assert y1_.max() / y1.max() == 10, "Scaling failed" + assert y2_.max() / y2.max() == 10, "Scaling failed" + + +def test_Lorentzian(): + post_sim = sp.SignalProcessor(operations=PS_1) + data = post_sim.apply_operations(data=sim.methods[0].simulation) + _, y0, y1, y2 = data.to_list() + + sigma = 200 + test = (sigma / 2) / (np.pi * (freqHz ** 2 + (sigma / 2) ** 2)) + + assert np.allclose(y1, y2) + assert np.all(y0 != y1) + assert np.allclose( + test / test.max(), y0 / y0.max(), atol=1e-04 + ), "Lorentzian apodization amplitude failed" + + +def test_Gaussian(): + post_sim = sp.SignalProcessor(operations=PS_2) + data = post_sim.apply_operations(data=sim.methods[0].simulation) + _, y0, y1, _ = data.to_list() + + sigma = 20 + test = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((freqHz / sigma) ** 2) / 2) + + assert np.allclose(y0, y1) + assert np.allclose( + test / test.max(), y0 / y0.max(), atol=1e-04 + ), "Gaussian apodization amplitude failed" + + # test None for dep_var_indx + post_sim = sp.SignalProcessor(operations=PS_3) + data = post_sim.apply_operations(data=sim.methods[0].simulation) + _, y0, y1, y2 = data.to_list() + + sigma = 20 + test = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((freqHz / sigma) ** 2) / 2) + + assert np.allclose(y0, y1) + assert np.allclose(y0, y2) + assert np.allclose( + test / test.max(), y0 / y0.max(), atol=1e-04 + ), "Gaussian apodization amplitude failed" + + +def test_scale_class(): + # direct initialization + a = sp.Scale(factor=200) + + assert a.factor == 200 + assert a.property_units == {} + + # class to dict with units + dict_ = a.to_dict_with_units() + + assert dict_ == { + "function": "Scale", + "factor": 200.0, + } + + # read from dictionary + b = sp.Scale.parse_dict_with_units(dict_) + + assert a == b + + +def test_Exponential_class(): + # direct initialization + a = apo.Exponential(FWHM=200, dim_indx=0, dep_var_indx=0) + + assert a.FWHM == 200 + assert a.property_units == {"FWHM": "Hz"} + assert a.dim_indx == 0 + assert a.dep_var_indx == 0 + + # class to dict with units + dict_ = a.to_dict_with_units() + + assert dict_ == { + "function": "apodization", + "type": "Exponential", + "FWHM": "200.0 Hz", + "dim_indx": 0, + "dep_var_indx": 0, + } + + # read from dictionary + b = apo.Exponential.parse_dict_with_units(dict_) + + assert a == b + + +def test_Gaussian_class(): + # direct initialization + a = apo.Gaussian(sigma=200, dim_indx=0, dep_var_indx=0) + + assert a.sigma == 200 + assert a.property_units == {"sigma": "Hz"} + assert a.dim_indx == 0 + assert a.dep_var_indx == 0 + + # class to dict with units + dict_ = a.to_dict_with_units() + + assert dict_ == { + "function": "apodization", + "type": "Gaussian", + "sigma": "200.0 Hz", + "dim_indx": 0, + "dep_var_indx": 0, + } + + # read from dictionary + b = apo.Gaussian.parse_dict_with_units(dict_) + + assert a == b diff --git a/src/mrsimulator/signal_processing/tests/test_signal_processing.py b/src/mrsimulator/signal_processing/tests/test_signal_processing.py new file mode 100644 index 000000000..9308fc31c --- /dev/null +++ b/src/mrsimulator/signal_processing/tests/test_signal_processing.py @@ -0,0 +1,44 @@ +# -*- coding: utf-8 -*- +import csdmpy as cp +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo +import numpy as np +import pytest + + +def test_01(): + post_sim = sp.SignalProcessor() + operations = [ + sp.IFFT(), + apo.Gaussian(sigma=12, dim_indx=0, dep_var_indx=0), + sp.FFT(), + ] + + post_sim.operations = operations + + with pytest.raises(ValueError, match="The data must be a CSDM object."): + post_sim.apply_operations([]) + + data = cp.as_csdm(np.arange(20)) + post_sim.apply_operations(data) + + # to dict with units + dict_ = post_sim.to_dict_with_units() + assert dict_ == { + "operations": [ + {"dim_indx": 0, "function": "IFFT"}, + { + "function": "apodization", + "type": "Gaussian", + "sigma": "12.0 Hz", + "dim_indx": 0, + "dep_var_indx": 0, + }, + {"dim_indx": 0, "function": "FFT"}, + ], + } + + # parse dict with units + post_sim_2 = sp.SignalProcessor.parse_dict_with_units(dict_) + + assert post_sim.operations == post_sim_2.operations diff --git a/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py new file mode 100644 index 000000000..a9e1fb74c --- /dev/null +++ b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py @@ -0,0 +1,104 @@ +# -*- coding: utf-8 -*- +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo +import mrsimulator.spectral_fitting as sf +from mrsimulator import Simulator +from mrsimulator import SpinSystem + + +def test_01(): + # str_to_html + str1 = "spin_systems[0].sites[0].isotropic_chemical_shift" + str1_encoded = "sys_0_site_0_isotropic_chemical_shift" + + str2 = "spin_systems[0].sites[0].quadrupolar.Cq" + str2_encoded = "sys_0_site_0_quadrupolar_Cq" + + assert sf._str_encode(str1) == str1_encoded + assert sf._str_encode(str2) == str2_encoded + + +def test_02(): + # html_to_str + str1 = "spin_systems[0].sites[0].isotropic_chemical_shift" + str1_encoded = "sys_0_site_0_isotropic_chemical_shift" + + str2 = "spin_systems[0].sites[0].quadrupolar.Cq" + str2_encoded = "sys_0_site_0_quadrupolar_Cq" + + assert sf._str_decode(str1_encoded) == str1 + assert sf._str_decode(str2_encoded) == str2 + + +def test_03(): + # post_sim_LMFIT_params + op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), + sp.Scale(factor=10), + ] + post_sim = sp.SignalProcessor(data=None, operations=op_list) + + params = sf._post_sim_LMFIT_params(post_sim) + + val = params.valuesdict() + + assert val["operation_1_Exponential"] == 100 + assert val["operation_3_Scale"] == 10 + + +def test_04(): + # update_post_sim_from_LMFIT_params + op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), + sp.Scale(factor=10), + ] + post_sim = sp.SignalProcessor(data=None, operations=op_list) + + params = sf._post_sim_LMFIT_params(post_sim) + + params["operation_1_Exponential"].value = 10 + params["operation_3_Scale"].value = 1 + + sf._update_post_sim_from_LMFIT_params(params, post_sim) + + assert post_sim.operations[1].FWHM == 10 + assert post_sim.operations[3].factor == 1 + + +def test_5(): + # make_LMFIT_parameters + sim = Simulator() + + H = { + "isotope": "1H", + "isotropic_chemical_shift": "10 ppm", + "shielding_symmetric": {"zeta": "5 ppm", "eta": 0.1}, + } + spin_system = {"sites": [H], "abundance": "100%"} + system_object = SpinSystem.parse_dict_with_units(spin_system) + sim.spin_systems += [system_object] + + op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), + sp.Scale(factor=10), + ] + post_sim = sp.SignalProcessor(data=None, operations=op_list) + + params = sf.make_LMFIT_parameters(sim, post_sim) + + valuesdict = { + "sys_0_site_0_isotropic_chemical_shift": 10, + "sys_0_site_0_shielding_symmetric_zeta": 5, + "sys_0_site_0_shielding_symmetric_eta": 0.1, + "sys_0_abundance": 100, + "operation_1_Exponential": 100, + "operation_3_Scale": 10, + } + + assert params.valuesdict() == valuesdict, "Parameter creation failed" diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 4dd4344e0..4a9d8b028 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -16,8 +16,6 @@ from .config import ConfigSimulator -# from mrsimulator.post_simulation import PostSimulator - __author__ = "Deepansh J. Srivastava" __email__ = "deepansh2012@gmail.com" @@ -57,7 +55,7 @@ class Simulator(BaseModel): >>> from mrsimulator.methods import BlochDecayCentralTransitionSpectrum >>> sim.methods = [ ... BlochDecaySpectrum(channels=['17O'], spectral_width=50000), - ... BlochDecayCentralTransitionSpectrum(channels=['17O'], spectral_width=50000) + ... BlochDecayCentralTransitionSpectrum(channels=['17O']) ... ] config: :ref:`config_api` object or equivalent dict object (optional). @@ -128,7 +126,6 @@ class Simulator(BaseModel): description: str = None spin_systems: List[SpinSystem] = [] methods: List[Method] = [] - # post_simulation: List[PostSimulator] = [] config: ConfigSimulator = ConfigSimulator() indexes = [] @@ -529,39 +526,3 @@ def _update_name_description_application(self, obj, index): "spin_systems": [self.spin_systems[index].to_dict_with_units()] } } - - def apodize(self, dimension=0, method=0, **kwargs): - if self.config.decompose_spectrum is False: - self.config.decompose_spectrum = True - self.run(method_index=method) - - csdm = self.methods[method].simulation - - for dim in csdm.dimensions: - dim.to("Hz", "nmr_frequency_ratio") - - for i, apodization_filter in enumerate(self.post_simulation): - apodization_filter.apodization[0]._apodize(csdm, i) - - for dim in csdm.dimensions: - dim.to("ppm", "nmr_frequency_ratio") - - # apodization_filter = Apodization( - # self.methods[method].simulation, dimension=dimension - # ) - # return apodization_filter.apodize(fn, **kwargs) - - # csdm = self.simulation - # for dim in csdm.dimensions: - # dim.to("Hz", "nmr_frequency_ratio") - # apo = self.post_simulation.apodization - - # sum_ = 0 - - # for item in apo: - # sum_ += item.apodize(csdm) - - # for dim in csdm.dimensions: - # dim.to("ppm", "nmr_frequency_ratio") - - # return self.post_simulation.scale * sum_ diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 6bcd44c2e..4228a4b3d 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- -import numpy as np +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo from mrsimulator import Simulator try: @@ -13,8 +14,32 @@ __author__ = "Maxwell C Venetos" __email__ = "maxvenetos@gmail.com" +START = "sys_" +ENCODING_PAIRS = [ + ["spin_systems[", START], + ["].sites[", "_site_"], + ["].isotropic_chemical_shift", "_isotropic_chemical_shift"], + ["].shielding_symmetric.", "_shielding_symmetric_"], + ["].quadrupolar.", "_quadrupolar_"], + ["].abundance", "_abundance"], + ["methods[", "METHODS_"], # why does methods need to be parameterized? + ["].post_simulation", "_POST_SIM_"], + [".scale", "scale"], + [".apodization[", "APODIZATION_"], + ["].args", "_args"], +] + +EXCLUDE = [ + "property_units", + "isotope", + "name", + "label", + "description", + "transition_pathways", +] -def _str_to_html(my_string): + +def _str_encode(my_string): """ LMFIT Parameters class does not allow for names to include special characters. This function converts '[', ']', and '.' to their HTML numbers to comply with @@ -25,26 +50,13 @@ def _str_to_html(my_string): Returns: String object. - """ - my_string = my_string.replace("spin_systems[", "ISO_") - my_string = my_string.replace("].sites[", "_SITES_") - my_string = my_string.replace( - "].isotropic_chemical_shift", "_isotropic_chemical_shift" - ) - my_string = my_string.replace("].shielding_symmetric.", "_shielding_symmetric_") - my_string = my_string.replace("].quadrupolar.", "_quadrupolar_") - my_string = my_string.replace("].abundance", "_abundance") - my_string = my_string.replace("methods[", "METHODS_") # - my_string = my_string.replace("].post_simulation", "_POST_SIM_") - my_string = my_string.replace(".scale", "scale") - my_string = my_string.replace(".apodization[", "APODIZATION_") - my_string = my_string.replace("].args", "_args") - + for item in ENCODING_PAIRS: + my_string = my_string.replace(*item) return my_string -def _html_to_string(my_string): +def _str_decode(my_string): """ Converts the HTML numbers to '[', ']', and '.' to allow for execution of the parameter name to update the simulator. @@ -54,23 +66,9 @@ def _html_to_string(my_string): Returns: String Object. - """ - my_string = my_string.replace("ISO_", "spin_systems[") - my_string = my_string.replace("_SITES_", "].sites[") - my_string = my_string.replace( - "_isotropic_chemical_shift", "].isotropic_chemical_shift" - ) - my_string = my_string.replace("_shielding_symmetric_", "].shielding_symmetric.") - my_string = my_string.replace("_quadrupolar_", "].quadrupolar.") - my_string = my_string.replace("_abundance", "].abundance") - - my_string = my_string.replace("METHODS_", "methods[") # - my_string = my_string.replace("_POST_SIM_", "].post_simulation") - my_string = my_string.replace("scale", ".scale") - my_string = my_string.replace("APODIZATION_", ".apodization[") - my_string = my_string.replace("_args", "].args") - + for item in ENCODING_PAIRS: + my_string = my_string.replace(*item[::-1]) return my_string @@ -84,21 +82,10 @@ def _list_of_dictionaries(my_list): Returns: List Object. - """ return [item.dict() for item in my_list] -exclude = [ - "property_units", - "isotope", - "name", - "label", - "description", - "transition_pathways", -] - - def _traverse_dictionaries(dictionary, parent="spin_systems"): """ Parses through the dictionary objects contained within the simulator object in @@ -106,50 +93,125 @@ def _traverse_dictionaries(dictionary, parent="spin_systems"): Args: dictionary: A dictionary or list object of the SpinSystem attributes from a - simulation object + simulation object parent: a string object used to create the addresses of the SpinSystem - attributes. + attributes. Returns: List Object. - """ name_list = [] if isinstance(dictionary, dict): for key, vals in dictionary.items(): - if key not in exclude and vals is not None: + if key not in EXCLUDE and vals is not None: if isinstance(vals, (dict, list)): name_list += _traverse_dictionaries( - vals, _str_to_html(f"{parent}.{key}") + vals, _str_encode(f"{parent}.{key}") ) else: - name_list += [_str_to_html(f"{parent}.{key}")] + name_list += [_str_encode(f"{parent}.{key}")] elif isinstance(dictionary, list): for i, items in enumerate(dictionary): - name_list += _traverse_dictionaries(items, _str_to_html(f"{parent}[{i}]")) - - # else: - # name_list += [_str_to_html(f"{parent}.{dictionary}")] + name_list += _traverse_dictionaries(items, _str_encode(f"{parent}[{i}]")) return name_list -def make_fitting_parameters(sim, exclude_key=None): +def _post_sim_LMFIT_params(post_sim): + """ + Creates an LMFIT Parameters object for SignalProcessor operations + involved in spectrum fitting + + Args: + post_sim: SignalProcessor object + + Returns: + Parameters object + """ + temp_dict = {} + # for item in post_sim.operations: + # prepend = f"DEP_VAR_{item.dependent_variable}_" + for i, operation in enumerate(post_sim.operations): + if isinstance(operation, apo.Gaussian): + identifier = f"operation_{i}_Gaussian" + arg = operation.sigma + temp_dict[f"{identifier}"] = arg + elif isinstance(operation, apo.Exponential): + identifier = f"operation_{i}_Exponential" + arg = operation.FWHM + temp_dict[f"{identifier}"] = arg + elif isinstance(operation, sp.Scale): + identifier = f"operation_{i}_Scale" + arg = operation.factor + temp_dict[f"{identifier}"] = arg + + params = Parameters() + for key, val in temp_dict.items(): + params.add(name=key, value=val) + + return params + + +def _update_post_sim_from_LMFIT_params(params, post_sim): + """Updates SignalProcessor operation arguments from an LMFIT Parameters object + + Args: + params: LMFIT Parameters object + post_sim: SignalProcessor object """ - Parses through the fitting parameter list to create LMFIT parameters used for - fitting. + temp_dict = {} + arg_dict = {"Gaussian": "sigma", "Exponential": "FWHM", "Scale": "factor"} + for param in params: + # iterating through the parameter list looking for only DEP_VAR + # (ie post_sim params) + if "operation_" in param: + # splitting parameter name to obtain + # Dependent variable index (var) + # index of operation in the operation list (opIndex) + # arg value for the operation (val) + split_name = param.split("_") + # var = split_name[split_name.index("VAR") + 1] + opIndex = split_name[split_name.index("operation") + 1] + val = params[param].value + # creating a dictionary of operations and arguments for each dependent + # variable + # if f"DepVar_{var}" not in temp_dict.keys(): + # temp_dict[f"DepVar_{var}"] = {} + temp_dict[f"{opIndex}_{split_name[-1]}"] = val + + # iterate through list of operation lists + # for item in post_sim.operations: + # iterating through dictionary with corresponding dependent variable index + for operation, val in temp_dict.items(): + # creating assignment strings to create the correct address for updating each + # operation + split = operation.split("_") + # dep_var_operation_list = f"post_sim.operations[{item.dependent_variable}]" + operation_val_update = f"post_sim.operations[{split[0]}].{arg_dict[split[-1]]}" + assignment = f"={val}" + exec(operation_val_update + assignment) + + +def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): + """ + Parses the Simulator and PostSimulator objects for a list of LMFIT parameters. + The parameter name is generated using the following syntax: + + ``sys_i_site_j_attribute1_attribute2`` + + for spin-system attribute with signature sys[i].sites[j].attribute1.attribute2 Args: sim: a Simulator object. + post_sim: a SignalProcessor object Returns: LMFIT Parameters object. - """ if not FOUND_LMFIT: error = ( f"The helper function {__name__} requires 'lmfit' module to create lmfit " - "paramters. Please install the lmfit module using\n'pip install lmfit'.", + r"paramters. Please install the lmfit module using\n'pip install lmfit'.", ) raise ImportError(error) @@ -158,110 +220,113 @@ def make_fitting_parameters(sim, exclude_key=None): params = Parameters() temp_list = _traverse_dictionaries(_list_of_dictionaries(sim.spin_systems)) - for i in range(len(sim.methods)): - if sim.methods[i].post_simulation is not None: - parent = f"methods[{i}].post_simulation" - - temp_list += [_str_to_html(parent + ".scale")] - # temp_list += [ - # item - # for item in _traverse_dictionaries( - # sim.methods[0].post_simulation, parent=parent - # ) - # if "scale" in item - # ] - if sim.methods[i].post_simulation.apodization is not None: - for j in range(len(sim.methods[i].post_simulation.apodization)): - temp_list.append(_str_to_html(parent + f".apodization[{j}].args")) + # get total abundance scaling factor length = len(sim.spin_systems) - abundance = 0 - last_abund = f"{length - 1}_abundance" - expression = "100" - for i in range(length - 1): - expression += f"-ISO_{i}_abundance" - for i in range(length): - abundance += eval("sim." + _html_to_string(f"spin_systems[{i}].abundance")) + abundance_scale = 100 / sum([sim.spin_systems[i].abundance for i in range(length)]) + # expression for the last abundance. + last_abund = f"{length - 1}_abundance" + expression = "-".join([f"{START}{i}_abundance" for i in range(length - 1)]) + expression = "100" if expression == "" else f"100-{expression}" for items in temp_list: - if "_eta" in items or "abundance" in items and last_abund not in items: - if "_eta" in items: - params.add( - name=items, - value=eval("sim." + _html_to_string(items)), - min=0, - max=1, - ) - if "abundance" in items: - params.add( - name=items, - value=eval("sim." + _html_to_string(items)) / abundance * 100, - min=0, - max=100, - ) + if "_eta" in items: + params.add( + name=items, value=eval("sim." + _str_decode(items)), min=0, max=1, + ) + # last_abund should come before abundance elif last_abund in items: params.add( name=items, - value=eval("sim." + _html_to_string(items)), + value=eval("sim." + _str_decode(items)), min=0, max=100, expr=expression, ) + elif "abundance" in items: + params.add( + name=items, + value=eval("sim." + _str_decode(items)) * abundance_scale, + min=0, + max=100, + ) else: - value = eval("sim." + _html_to_string(items)) - if type(value) == list: - params.add(name=items, value=value[0]) - else: - params.add(name=items, value=value) + value = eval("sim." + _str_decode(items)) + params.add(name=items, value=value) + + if post_sim is None: + return params + + if not isinstance(post_sim, sp.SignalProcessor): + raise ValueError( + f"Expecting a `SignalProcessor` object, found {type(post_sim).__name__}." + ) + temp_params = _post_sim_LMFIT_params(post_sim) + for item in temp_params: + params.add(name=item, value=temp_params[item].value) + # params.add_many(temp_params) return params -def min_function(params, sim, apodization_function=None): +def LMFIT_min_function(params, sim, post_sim=None): """ - The simulation routine to establish how the parameters will update the simulation. + The simulation routine to calculate the vector difference between simulation and + experiment based on the parameters update. Args: params: Parameters object containing parameters to vary during minimization. - data: a CSDM object of the data to fit the simulation to. - sim: Simulator object to be fit to data. Initialized with arbitrary fitting - parameters. - apodization_function: A string indicating the apodization function to use. - Currently "Gaussian" and "Lorentzian" are supported. + sim: Simulator object used in the simulation. Initialized with guess fitting + parameters. + post_sim: PostSimulator object used in the simulation. Initialized with guess + fitting parameters. Returns: Array of the differences between the simulation and the experimental data. - """ if not isinstance(params, Parameters): raise ValueError( f"Expecting a `Parameters` object, found {type(params).__name__}." ) - # if not isinstance(data, cp.CSDM): - # raise ValueError(f"Expecting a `CSDM` object, found {type(data).__name__}.") + if not isinstance(post_sim, sp.SignalProcessor) or post_sim is None: + raise ValueError( + f"Expecting a `SignalProcessor` object, found {type(post_sim).__name__}." + ) + if not isinstance(sim, Simulator): raise ValueError(f"Expecting a `Simulator` object, found {type(sim).__name__}.") - # intensity_data = data.dependent_variables[0].components[0].real values = params.valuesdict() for items in values: - if "args" not in items: - nameString = "sim." + _html_to_string(items) + if "operation_" not in items: + nameString = "sim." + _str_decode(items) executable = f"{nameString} = {values[items]}" exec(executable) - else: - nameString = "sim." + _html_to_string(items) - executable = f"{nameString} = [{values[items]}]" - exec(executable) + elif "operation_" in items and post_sim is not None: + _update_post_sim_from_LMFIT_params(params, post_sim) sim.run() - residual = np.asarray([]) + processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) + # residual = np.asarray([]) + + if sim.config.decompose_spectrum == "spin_system": + datum = 0 + for decomposed_datum in processed_data.dependent_variables: + datum += decomposed_datum.components[0] + # datum = [sum(i) for i in zip(datum, decomposed_datum)] + else: + datum = processed_data.dependent_variables[0].components[0] + + return ( + sim.methods[0].experiment.dependent_variables[0].components[0].real - datum.real + ) - for i, method in enumerate(sim.methods): - y_factored = method.apodize().real - residual = np.append( - residual, - method.experiment.dependent_variables[0].components[0].real - y_factored, - ) + # MULTIPLE EXPERIMENTS + # for i, method in enumerate(sim.methods): + # y_factored = method.apodize().real + # residual = np.append( + # residual, + # method.experiment.dependent_variables[0].components[0].real - y_factored, + # ) - return residual + # return residual diff --git a/src/mrsimulator/spin_system/__init__.py b/src/mrsimulator/spin_system/__init__.py index 031700285..d927d76d0 100644 --- a/src/mrsimulator/spin_system/__init__.py +++ b/src/mrsimulator/spin_system/__init__.py @@ -121,7 +121,8 @@ class SpinSystem(Parseable): to the computation. .. seealso:: - `Fitting example <./../auto_examples/Fitting/plot_2_mrsimFitExample_O17.html>`_ + `Fitting example + <./../auto_examples/Fitting/plot_2_mrsimFitExample_O17.html>`_ """ name: str = None diff --git a/src/mrsimulator/spin_system/tensors.py b/src/mrsimulator/spin_system/tensors.py index 8feabac36..eb8097da9 100644 --- a/src/mrsimulator/spin_system/tensors.py +++ b/src/mrsimulator/spin_system/tensors.py @@ -79,7 +79,7 @@ class SymmetricTensor(Parseable): Example ------- - >>> shielding = SymmetricTensor(zeta=10, eta=0.1, alpha=0.15, beta=3.1415, gamma=2.1) + >>> shielding = SymmetricTensor(zeta=10, eta=0.1, alpha=0.15, beta=3.14, gamma=2.1) >>> efg = SymmetricTensor(Cq=10e6, eta=0.5, alpha=1.5, beta=1.1451, gamma=0) """ diff --git a/tests/test_apodization.py b/tests/test_apodization.py deleted file mode 100644 index 48d9bacb7..000000000 --- a/tests/test_apodization.py +++ /dev/null @@ -1,65 +0,0 @@ -# -*- coding: utf-8 -*- -"""Apodization test""" -import numpy as np -from mrsimulator import Simulator -from mrsimulator import SpinSystem -from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.post_simulation import PostSimulator - - -sim = Simulator() -the_site = {"isotope": "1H", "isotropic_chemical_shift": "0 ppm"} -the_spin_system = {"name": "site A", "sites": [the_site], "abundance": "80%"} -spin_system_object = SpinSystem.parse_dict_with_units(the_spin_system) -sim.spin_systems += [spin_system_object] - -method_1 = BlochDecaySpectrum( - channels=["1H"], - magnetic_flux_density=9.4, - rotor_angle=0, - rotor_frequency=0, - spectral_dimensions=[ - {"count": 4096, "spectral_width": 25000, "reference_offset": 0} - ], -) - - -PS_1 = PostSimulator( - scale=1, apodization=[{"args": [200], "function": "Lorentzian", "dimension": 0}] -) - -PS_2 = PostSimulator( - scale=1, apodization=[{"args": [20], "function": "Gaussian", "dimension": 0}] -) - -sim.methods += [method_1, method_1] -sim.run() - - -freqHz = sim.methods[0].spectral_dimensions[0].coordinates_Hz() - - -def test_Lorentzian(): - sim.methods[0].post_simulation = PS_1 - - sigma = 200 - test = (sigma / 2) / (np.pi * (freqHz ** 2 + (sigma / 2) ** 2)) - - x, y = sim.methods[0].simulation.to_list() - y = sim.methods[0].apodize().real - - assert np.allclose( - test / test.max(), y / y.max(), atol=1e-04 - ), "Lorentzian apodization amplitude failed" - - -def test_Gaussian(): - sim.methods[0].post_simulation = PS_2 - sigma = 20 - test = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((freqHz / sigma) ** 2) / 2) - x, y = sim.methods[1].simulation.to_list() - y = sim.methods[0].apodize().real - - assert np.allclose( - test / test.max(), y / y.max(), atol=1e-04 - ), "Gaussian apodization amplitude failed" diff --git a/tests/test_parameters.py b/tests/test_parameters.py deleted file mode 100644 index b041f9871..000000000 --- a/tests/test_parameters.py +++ /dev/null @@ -1,42 +0,0 @@ -# -*- coding: utf-8 -*- -"""Parameter test""" -from mrsimulator import Simulator -from mrsimulator import Site -from mrsimulator import SpinSystem -from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.spectral_fitting import make_fitting_parameters -from mrsimulator.spin_system.tensors import SymmetricTensor as st - - -sim = Simulator() - -H = Site( - isotope="1H", isotropic_chemical_shift=10, shielding_symmetric=st(zeta=5, eta=0.1) -) -spin_system = SpinSystem(name="H1", sites=[H], abundance=100) - -method = BlochDecaySpectrum( - channels=["1H"], - magnetic_flux_density=9.4, - rotor_frequency=14000, - spectral_dimensions=[ - {"count": 2048, "spectral_width": 20000, "reference_offset": 0} - ], -) - - -sim.spin_systems += [spin_system] -sim.methods += [method] - -params = make_fitting_parameters(sim) - -valuesdict = { - "ISO_0_SITES_0_isotropic_chemical_shift": 10, - "ISO_0_SITES_0_shielding_symmetric_zeta": 5, - "ISO_0_SITES_0_shielding_symmetric_eta": 0.1, - "ISO_0_abundance": 100, -} - - -def test_param(): - assert params.valuesdict() == valuesdict, "Parameter creation failed"