From 25a01c051aa90f8f24bae96eb703e3d93e903df1 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Mon, 15 Jun 2020 16:41:17 -0400 Subject: [PATCH 01/26] post simulation signal processing --- .flake8 | 2 +- src/mrsimulator/method/__init__.py | 30 ---- src/mrsimulator/post_simulation.py | 225 +++++++++++++++++++------- src/mrsimulator/simulator/__init__.py | 37 ----- tests/test_apodization.py | 53 ++++-- 5 files changed, 205 insertions(+), 142 deletions(-) diff --git a/.flake8 b/.flake8 index 969b66bdb..e598e9d8e 100644 --- a/.flake8 +++ b/.flake8 @@ -1,5 +1,5 @@ [flake8] -ignore = E402 C901 E501 E265 +ignore = E402 #C901 E265 E402 E501 exclude = .eggs, *.egg,build, src/mrsimulator/__init__.py filename = *.pyx, *py max-line-length = 88 diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 05d379627..c83624076 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 @@ -55,14 +54,12 @@ class Method(Parseable): or a `ndarray `_ object holding the experimental measurement for the given method, if available. The default value is None. - 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 @@ -159,8 +156,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: @@ -282,28 +277,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 index bffe02242..e132cb58c 100644 --- a/src/mrsimulator/post_simulation.py +++ b/src/mrsimulator/post_simulation.py @@ -1,98 +1,207 @@ # -*- coding: utf-8 -*- """The Event class.""" +from typing import ClassVar +from typing import Dict from typing import List -from typing import Optional +from typing import Union +import csdmpy as cp import numpy as np from mrsimulator.util.parseable import Parseable -from numpy.fft import fft -from numpy.fft import fftshift -from numpy.fft import ifft -from numpy.fft import ifftshift +from pydantic import BaseModel -class Apodization(Parseable): - """The Apodization class.""" +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.__class__.__name__ + return my_dict + + @classmethod + def parse_dict_with_units(cls, py_dict): + 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) + + def operate(self, data): + return data + + +class Scale(abstractOperation): + """ + abstractOperation-based class for applying a scaling + factor to a dependent variable of simulation data + """ + + factor: float = 1 + + def operate(self, data, dep_var): + """ + 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 + """ + data_copy = data.copy() + factored = data_copy.dependent_variables[dep_var].components[0] * self.factor + data_copy.dependent_variables[dep_var].components[0] = factored + return data_copy + + +class Gaussian(abstractOperation): + """ + abstractOperation-based class for applying a Gaussian function + to a dependent variable of simulation data + + dimension: int. Data dimension to apply the function along + sigma: float. Standard deviation of Gaussian function + """ - function: str dimension: int = 0 - args: List[float] = [0] - fraction: Optional[float] = 1 + sigma: float = 0 - class Config: - validate_assignment = True + property_unit_types: ClassVar = {"sigma": ["time", "frequency"]} + property_default_units: ClassVar = {"sigma": ["s", "Hz"]} + property_units: Dict = {"sigma": ["s", "Hz"]} + + def operate(self, data, dep_var): + """ + 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 + """ + # unit = self.property_units["sigma"] + data_copy = data.copy() + # data_copy.dimensions[self.dimension].coordinates.to(unit).value + x = data_copy.dimensions[self.dimension].coordinates.value + data_copy.dependent_variables[dep_var].components[0] *= np.exp( + -((x * self.sigma * np.pi) ** 2) * 2 + ) + return data_copy - def apodize(self, csdm, **kwargs): - """Returns the result of passing the selected apodization function . - Args: - csdm: simulation object +class Exponential(abstractOperation): + """ + abstractOperation-based class for applying an exponential + Lorentzian function to a dependent variable of simulation data + + dimension: int. Data dimension to apply the function along + Lambda: float. Width parameter + """ - Returns: - A Numpy array + dimension: int = 0 + Lambda: float = 0 + + property_unit_types: ClassVar = {"Lambda": ["time", "frequency"]} + property_default_units: ClassVar = {"Lambda": ["s", "Hz"]} + property_units: Dict = {"Lambda": ["s", "Hz"]} + + def operate(self, data, dep_var): """ - function_mapping = {"Gaussian": Gaussian, "Lorentzian": Lorentzian} + Applies the operation for which the class is named for. - 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 + data: CSDM object + dep_var: int. The index of the dependent variable to apply operation to + """ + # unit = self.property_units["Lambda"] + data_copy = data.copy() + # x = data.dimensions[self.dimension].coordinates.to(unit).value + x = data.dimensions[self.dimension].coordinates.value + data_copy.dependent_variables[dep_var].components[0] *= np.exp( + -self.Lambda * np.pi * np.abs(x) ) + return data_copy - 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) +class IFFT(abstractOperation): + """ + abstractOperation-based class for applying an inverse + Fourier transform to a dependent variable of simulation data - appodized = TimeDomain * fapp - fft_output = fftshift(fft(appodized, axis=axis), axes=axis).real + dimension: int. Data dimension to apply the function along - # csdm.dependent_variables[index].components[0] = ( - # self.fraction * phase.conj() * fft_output - # ) + """ - return self.fraction * phase.conj() * fft_output + dimension: int = 0 + def operate(self, data, dep_var): + """ + Applies the operation for which the class is named for. -def Lorentzian(csdm, axis, arg, **kwargs): - """Lorentzian apodization function. + data: CSDM object + dep_var: int. The index of the dependent variable to apply operation to + """ + data_copy = data.copy() + return data_copy.fft(axis=self.dimension) - Args: - Self: simulation object - arg: list with one entry. The full-width-half-max in Hz of the Lorentzian function - Returns: - A Numpy array +class FFT(abstractOperation): """ - inv_x = csdm.dimensions[axis].reciprocal_coordinates().to("s").value - return np.exp(-arg[0] * np.pi * np.abs(inv_x)) + abstractOperation-based class for applying a + Fourier transform to a dependent variable of simulation data + dimension: int. Data dimension to apply the function along -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 + dimension: int = 0 - Returns: - A Numpy array + def operate(self, data, dep_var): + """ + 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 + """ + data_copy = data.copy() + return data_copy.fft(axis=self.dimension) + + +class OperationList(BaseModel): + """ + Container class to hold a list of operations to be applied to a given + dependent variable in the data + + dependent_variable: The dependent variable of the data to apply operations to + operations: List of abstractOperation based operations to sequentially apply to data. """ - inv_x = csdm.dimensions[axis].reciprocal_coordinates().to("s").value - return np.exp(-((inv_x * arg[0] * np.pi) ** 2) * 2) + dependent_variable: int = 0 + operations: List[abstractOperation] -class PostSimulator(Parseable): - """scale: scaling factor + +class SignalProcessor(BaseModel): + """ + Signal processing class to apply lists of various operations to individual dependent variables + of the data. + + data: CSDM object. From simulation + operations: List of operation lists """ - scale: Optional[float] = 1.0 - apodization: List[Apodization] = [] - truncation_artifact: Optional[bool] = False - dead_time: Optional[float] = 0.0 + data: Union[cp.CSDM, np.ndarray] = None + operations: List[OperationList] class Config: validate_assignment = True arbitrary_types_allowed = True + + def apply_operations(self, **kwargs): + op_list = self.operations + + for item in op_list: + dep_var = item.dependent_variable + for filters in item.operations: + self.data = filters.operate(self.data, dep_var=dep_var) + + return self.data diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 1fd8a2729..8b7b4fd8e 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -58,7 +58,6 @@ def __eq__(self, other): self.description == other.description, self.spin_systems == other.spin_systems, self.methods == other.methods, - # self.post_simulation == other.post_simulation, self.config == other.config, ] if np.all(check): @@ -394,39 +393,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/tests/test_apodization.py b/tests/test_apodization.py index 740dbe7ba..3ee3743b2 100644 --- a/tests/test_apodization.py +++ b/tests/test_apodization.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """Apodization test""" +import mrsimulator.post_simulation as ps import numpy as np from mrsimulator import Simulator from mrsimulator import SpinSystem from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.post_simulation import PostSimulator - +from mrsimulator.post_simulation import SignalProcessor sim = Simulator() the_site = {"isotope": "1H", "isotropic_chemical_shift": "0 ppm"} @@ -23,42 +23,63 @@ ], ) +PS_0 = {"dependent_variable": 0, "operations": [ps.Scale(factor=10)]} -PS_1 = PostSimulator( - scale=1, apodization=[{"args": [200], "function": "Lorentzian", "dimension": 0}] -) - -PS_2 = PostSimulator( - scale=1, apodization=[{"args": [20], "function": "Gaussian", "dimension": 0}] -) +PS_1 = { + "dependent_variable": 0, + "operations": [ + ps.IFFT(dimension=0), + ps.Exponential(Lambda=200, dimension=0), + ps.FFT(dimension=0), + ], +} + +PS_2 = { + "dependent_variable": 0, + "operations": [ + ps.IFFT(dimension=0), + ps.Gaussian(sigma=20, dimension=0), + ps.FFT(dimension=0), + ], +} -sim.methods += [method_1, method_1] +sim.methods += [method_1] sim.run() freqHz = sim.methods[0].spectral_dimensions[0].coordinates_Hz +def test_scale(): + post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_0]) + post_sim.apply_operations() + x0, y0 = sim.methods[0].simulation.to_list() + x, y = post_sim.data.to_list() + + assert np.allclose(sum(y0 / y), 10), "Scaling failed" + + def test_Lorentzian(): sim.methods[0].post_simulation = PS_1 + post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_1]) + post_sim.apply_operations() + x, y = post_sim.data.to_list() 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 + post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_2]) + post_sim.apply_operations() + x, y = post_sim.data.to_list() + 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 From d8a80a9870c95f6be5fdf3cfe9efaf0e3d08961e Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Mon, 15 Jun 2020 16:44:56 -0400 Subject: [PATCH 02/26] minor text changes after merging --- src/mrsimulator/method/__init__.py | 1 - src/mrsimulator/simulator/__init__.py | 3 --- 2 files changed, 4 deletions(-) diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 364a81145..4e1948e16 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -102,7 +102,6 @@ class Method(Parseable): >>> bloch.description 'Huh!' - post_simulation: An optional dict with post-simulation parameters. """ name: str = None label: str = None diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index ef90c519c..92cd08853 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" @@ -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 = [] From bc52a81ef9ade4f56029f26c1300195637f9a069 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Tue, 16 Jun 2020 13:24:26 -0400 Subject: [PATCH 03/26] Test errors and unused features in spectral fit --- src/mrsimulator/post_simulation.py | 3 +-- src/mrsimulator/spectral_fitting.py | 30 ++++++++++++++--------------- tests/test_apodization.py | 3 +-- 3 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/mrsimulator/post_simulation.py b/src/mrsimulator/post_simulation.py index e132cb58c..5c10d6651 100644 --- a/src/mrsimulator/post_simulation.py +++ b/src/mrsimulator/post_simulation.py @@ -52,8 +52,7 @@ def operate(self, data, dep_var): dep_var: int. The index of the dependent variable to apply operation to """ data_copy = data.copy() - factored = data_copy.dependent_variables[dep_var].components[0] * self.factor - data_copy.dependent_variables[dep_var].components[0] = factored + data_copy.dependent_variables[dep_var].components[0] *= self.factor return data_copy diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 6bcd44c2e..708e9bb1f 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -158,21 +158,21 @@ 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")) + # 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")) length = len(sim.spin_systems) abundance = 0 diff --git a/tests/test_apodization.py b/tests/test_apodization.py index ebcd53e0b..a49f88de6 100644 --- a/tests/test_apodization.py +++ b/tests/test_apodization.py @@ -56,11 +56,10 @@ def test_scale(): x0, y0 = sim.methods[0].simulation.to_list() x, y = post_sim.data.to_list() - assert np.allclose(sum(y0 / y), 10), "Scaling failed" + assert y.max() / y0.max() == 10, "Scaling failed" def test_Lorentzian(): - sim.methods[0].post_simulation = PS_1 post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_1]) post_sim.apply_operations() x, y = post_sim.data.to_list() From adaf7656a4128621cadb182ebce03235986105a5 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Tue, 16 Jun 2020 14:04:20 -0400 Subject: [PATCH 04/26] fixed Method doctest --- src/mrsimulator/method/__init__.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 4e1948e16..87804d268 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -41,8 +41,8 @@ class Method(Parseable): Example ------- - >>> bloch = Method() - >>> bloch.channels = ['1H'] + >>> bloch = Method() + >>> bloch.channels = ['1H'] spectral_dimensions: list of :ref:`spectral_dim_api` or dict objects (optional). The number of spectral dimensions depends on the given method. For example, a @@ -106,8 +106,8 @@ class Method(Parseable): name: str = None label: str = None description: str = None - channels: List[str] - spectral_dimensions: List[SpectralDimension] + channels: List[str] = [] + spectral_dimensions: List[SpectralDimension] = [{}] simulation: Union[cp.CSDM, np.ndarray] = None experiment: Union[cp.CSDM, np.ndarray] = None From 304903f71f3843b66b11ca785f338904bb203774 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Tue, 16 Jun 2020 16:11:14 -0400 Subject: [PATCH 05/26] post_simulation data_copy efficiency --- src/mrsimulator/post_simulation.py | 47 +++++++++++++++++------------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/src/mrsimulator/post_simulation.py b/src/mrsimulator/post_simulation.py index 5c10d6651..d87160eae 100644 --- a/src/mrsimulator/post_simulation.py +++ b/src/mrsimulator/post_simulation.py @@ -38,7 +38,7 @@ def operate(self, data): class Scale(abstractOperation): """ - abstractOperation-based class for applying a scaling + Class for applying a scaling factor to a dependent variable of simulation data """ @@ -48,19 +48,22 @@ def operate(self, data, dep_var): """ Applies the operation for which the class is named for. + $f(\vec(x)) = scale*\vec(x)$ + data: CSDM object dep_var: int. The index of the dependent variable to apply operation to """ - data_copy = data.copy() - data_copy.dependent_variables[dep_var].components[0] *= self.factor - return data_copy + data.dependent_variables[dep_var].components[0] *= self.factor + return data class Gaussian(abstractOperation): """ - abstractOperation-based class for applying a Gaussian function + Class for applying a Gaussian function to a dependent variable of simulation data + $f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2}$ + dimension: int. Data dimension to apply the function along sigma: float. Standard deviation of Gaussian function """ @@ -80,20 +83,21 @@ def operate(self, data, dep_var): dep_var: int. The index of the dependent variable to apply operation to """ # unit = self.property_units["sigma"] - data_copy = data.copy() # data_copy.dimensions[self.dimension].coordinates.to(unit).value - x = data_copy.dimensions[self.dimension].coordinates.value - data_copy.dependent_variables[dep_var].components[0] *= np.exp( + x = data.dimensions[self.dimension].coordinates.value + data.dependent_variables[dep_var].components[0] *= np.exp( -((x * self.sigma * np.pi) ** 2) * 2 ) - return data_copy + return data class Exponential(abstractOperation): """ - abstractOperation-based class for applying an exponential + Class for applying an exponential Lorentzian function to a dependent variable of simulation data + $f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)}$ + dimension: int. Data dimension to apply the function along Lambda: float. Width parameter """ @@ -113,18 +117,17 @@ def operate(self, data, dep_var): dep_var: int. The index of the dependent variable to apply operation to """ # unit = self.property_units["Lambda"] - data_copy = data.copy() # x = data.dimensions[self.dimension].coordinates.to(unit).value x = data.dimensions[self.dimension].coordinates.value - data_copy.dependent_variables[dep_var].components[0] *= np.exp( + data.dependent_variables[dep_var].components[0] *= np.exp( -self.Lambda * np.pi * np.abs(x) ) - return data_copy + return data class IFFT(abstractOperation): """ - abstractOperation-based class for applying an inverse + Class for applying an inverse Fourier transform to a dependent variable of simulation data dimension: int. Data dimension to apply the function along @@ -140,13 +143,12 @@ def operate(self, data, dep_var): data: CSDM object dep_var: int. The index of the dependent variable to apply operation to """ - data_copy = data.copy() - return data_copy.fft(axis=self.dimension) + return data.fft(axis=self.dimension) class FFT(abstractOperation): """ - abstractOperation-based class for applying a + Class for applying a Fourier transform to a dependent variable of simulation data dimension: int. Data dimension to apply the function along @@ -162,8 +164,7 @@ def operate(self, data, dep_var): data: CSDM object dep_var: int. The index of the dependent variable to apply operation to """ - data_copy = data.copy() - return data_copy.fft(axis=self.dimension) + return data.fft(axis=self.dimension) class OperationList(BaseModel): @@ -198,9 +199,13 @@ class Config: def apply_operations(self, **kwargs): op_list = self.operations + copy_data = self.data.copy() + for item in op_list: dep_var = item.dependent_variable for filters in item.operations: - self.data = filters.operate(self.data, dep_var=dep_var) + copy_data = filters.operate(copy_data, dep_var=dep_var) + + self.data = copy_data - return self.data + return copy_data From 21d07afbce75277144b65ebade303b2dd0a47220 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Thu, 18 Jun 2020 15:55:05 -0400 Subject: [PATCH 06/26] fitting updates and fitting docs --- .../plot_1_mrsimFitExample_cuspidine.py | 141 +++++++---- .../Fitting/plot_2_mrsimFitExample_O17.py | 153 ++++++------ src/mrsimulator/post_simulation.py | 12 +- src/mrsimulator/spectral_fitting.py | 227 ++++++++++++------ 4 files changed, 337 insertions(+), 196 deletions(-) diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 67a86aee0..d21a9f502 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -6,7 +6,7 @@ .. 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 @@ -20,13 +20,14 @@ # The :math:`^{29}\text{Si}` tensor parameters were obtained from Hansen `et. al.` [#f1]_ # # 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 = {"weight": "light", "size": 9} +mpl.rc("font", **font) +mpl.rcParams["figure.figsize"] = [4.25, 3.0] -#%% +# %% # Next we will import `csdmpy `_ and loading the data file. import csdmpy as cp @@ -34,63 +35,104 @@ synthetic_experiment = cp.load(filename) synthetic_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") -x1, y1 = synthetic_experiment.to_list() - -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) +ax = plt.subplot(projection="csdm") +ax.plot(synthetic_experiment, color="black", linewidth=1) +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 SpinSystem +from mrsimulator import Simulator +# %% +# **Step 1** Create the site. -from mrsimulator import Simulator, SpinSystem, Site -from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.post_simulation import PostSimulator +site = { + "isotope": "29Si", + "isotropic_chemical_shift": "-80.0 ppm", + "shielding_symmetric": {"zeta": "-60 ppm", "eta": 0.6}, +} -S29 = Site( - isotope="29Si", - isotropic_chemical_shift=-80.0, - shielding_symmetric={"zeta": -60, "eta": 0.6}, -) +# %% +# **Step 2** Create the spin system for the site. + +spin_system = { + "name": "Si Site", + "description": "A 29Si site in cuspidine", + "sites": [site], # from the above code + "abundance": "100%", +} +system_object = SpinSystem.parse_dict_with_units(spin_system) + +# %% +# **Step 3** Create the Bloch Decay method. + +from mrsimulator.methods import BlochDecaySpectrum 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": 2046, "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}] -) +# %% +# The above method is set up to record the :math:`^{29}\text{Si}` resonances at the +# magic angle, spinning at 780 Hz and 7.1 T external magnetic flux density. +# The resonances are recorded over 25 kHz spectral width ofset by -5000 Hz and +# using 2046 points. +# %% +# **Step 4** 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.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 + +import mrsimulator.post_simulation as ps +from mrsimulator.post_simulation import SignalProcessor + +op_list = { + "dependent_variable": 0, + "operations": [ + ps.IFFT(dimension=0), + ps.Exponential(Lambda=200, dimension=0), + ps.FFT(dimension=0), + ps.Scale(factor=1), + ], +} + +post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) + + +# %% +# ** Step 7** Process and plot the spectrum. + +post_sim.apply_operations() -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) +ax = plt.subplot(projection="csdm") +ax.plot(post_sim.data, color="black", linewidth=1) +ax.invert_xaxis() plt.tight_layout() plt.show() @@ -105,15 +147,15 @@ #%% -from lmfit import Minimizer, Parameters +from lmfit import Minimizer, Parameters, fit_report 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.add(name="Lambda", value=200) +params.add(name="factor", value=1) params #%% @@ -132,31 +174,33 @@ def test_function(params, data, sim): 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 sim.run() - # here we apodize the signal to simulate line broadening - y = sim.methods[0].apodize().real - y_factored = y * values["scale"] + post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) + post_sim.operations[0].operations[3].factor = values["factor"] + post_sim.operations[0].operations[1].Lambda = values["Lambda"] + # here we apodize the signal to simulate line broadening + post_sim.apply_operations() - return intensity - y_factored + return intensity - post_sim.data.dependent_variables[0].components[0].real #%% # 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 magic angle spinning fitting is we must use a discrete minimization method +# such as 'powell' as the chemical shift varies discretely -#%% +# %% **Step 8** Perform minimization. minner = Minimizer(test_function, params, fcn_args=(synthetic_experiment, sim)) result = minner.minimize(method="powell") -result +print(fit_report(result)) #%% -# After the fit, we can plot the new simulated spectrum. +# **Step 9** Plot the fitted spectrum. plt.figsize = (4, 3) @@ -166,7 +210,6 @@ def test_function(params, data, sim): plt.plot(*(synthetic_experiment - residual).to_list(), "r", alpha=0.5, label="Fit") plt.plot(*residual.to_list(), alpha=0.5, label="Residual") -plt.xlim(x1.value.max(), x1.value.min()) plt.xlabel("Frequency / Hz") plt.grid(which="major", axis="both", linestyle="--") plt.legend() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 1ddb2e563..a4d32b4e8 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -19,11 +19,12 @@ # We use the :math:`^{17}\text{O}` tensor information from Grandinetti `et. al.` [#f5]_ # # 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 = {"weight": "light", "size": 9} +mpl.rc("font", **font) +mpl.rcParams["figure.figsize"] = [4.25, 3.0] #%% # Next we will import `csdmpy `_ and loading the data file. @@ -34,12 +35,9 @@ oxygen_experiment = cp.load(filename).real oxygen_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") -x1, y1 = oxygen_experiment.to_list() - -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) +ax = plt.subplot(projection="csdm") +ax.plot(oxygen_experiment, color="black", linewidth=1) +ax.invert_xaxis() plt.tight_layout() plt.show() @@ -49,24 +47,36 @@ # 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 +from mrsimulator import SpinSystem +from mrsimulator import Simulator -sim = Simulator() -O17_1 = Site( - isotope="17O", - isotropic_chemical_shift=60.0, - quadrupolar={"Cq": 4200000, "eta": 0.5}, -) -O17_2 = Site( - isotope="17O", isotropic_chemical_shift=40, quadrupolar={"Cq": 2400000, "eta": 0.18} -) +# %% +# **Step 1** Create the site. + +O17_1 = { + "isotope": "17O", + "isotropic_chemical_shift": "60.0 ppm", + "quadrupolar": {"Cq": "4200000 Hz", "eta": 0.5}, +} + +O17_2 = { + "isotope": "17O", + "isotropic_chemical_shift": "40.0 ppm", + "quadrupolar": {"Cq": "2400000 Hz", "eta": 0.5}, +} -spin_systems = [SpinSystem(sites=[site]) for site in [O17_1, O17_2]] +# %% +# **Step 2** Create the spin system for the site. +spin_system = {"sites": [O17_1, O17_2], "abundance": "100%"} # from the above code + +system_object = SpinSystem.parse_dict_with_units(spin_system) + +# %% +# **Step 3** Create the Bloch Decay method. + +from mrsimulator.methods import BlochDecayCentralTransitionSpectrum count = oxygen_experiment.dimensions[0].count increment = oxygen_experiment.dimensions[0].increment.to("Hz").value @@ -74,42 +84,62 @@ 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 } ], ) +# %% +# **Step 4** Create the Simulator object and add the method and spin-system objects. -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 = Simulator() +sim.spin_systems += [system_object] sim.methods += [method] -sim.methods[0].post_simulation = PS + sim.methods[0].experiment = oxygen_experiment -# 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 we will save the relevant transition pathways iso.transition_pathways = method.get_transition_pathways(iso).tolist() sim.run() -sim.methods[0].simulation.dimensions[0].to("ppm", "nmr_frequency_ratio") +# %% +# **Step 6** Create a SignalProcessor -x, y = sim.methods[0].simulation.to_list() -y = sim.methods[0].apodize().real +import mrsimulator.post_simulation as ps +from mrsimulator.post_simulation import SignalProcessor -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) +factor = oxygen_experiment.dependent_variables[0].components[0].max().real / 4 + +op_list = { + "dependent_variable": 0, + "operations": [ + ps.IFFT(dimension=0), + ps.Exponential(Lambda=100, dimension=0), + ps.FFT(dimension=0), + ps.Scale(factor=factor), + ], +} + +post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) + +# %% +# ** Step 7** Process and plot the spectrum. + +post_sim.apply_operations() + +ax = plt.subplot(projection="csdm") +ax.plot(post_sim.data, color="black", linewidth=1) +ax.invert_xaxis() plt.tight_layout() plt.show() @@ -121,40 +151,26 @@ # attributes may vary. To simplify the parameter list creation we will use the # :func:`~mrsimulator.spectral_fitting.make_fitting_parameters` -#%% +# %% +# **Step 8** Create a list of parameters to vary during fitting. 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. -# -# 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. - - -#%% +params = make_fitting_parameters(sim, post_sim) +params.pretty_print() +# %% +# **Step 9** Perform minimization. from mrsimulator.spectral_fitting import min_function -from lmfit import Minimizer +from lmfit import Minimizer, report_fit -minner = Minimizer(min_function, params, fcn_args=(sim, "Lorentzian")) +minner = Minimizer(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** Plot the fitted spectrum. plt.figsize = (4, 3) @@ -164,7 +180,6 @@ plt.plot(*(oxygen_experiment - residual).to_list(), "r", alpha=0.5, label="Fit") plt.plot(*residual.to_list(), alpha=0.5, label="Residual") -plt.xlim(x.value.max(), x.value.min()) plt.xlabel("$^{17}$O frequency / ppm") plt.grid(which="major", axis="both", linestyle="--") plt.legend() diff --git a/src/mrsimulator/post_simulation.py b/src/mrsimulator/post_simulation.py index d87160eae..ad153b20d 100644 --- a/src/mrsimulator/post_simulation.py +++ b/src/mrsimulator/post_simulation.py @@ -32,9 +32,6 @@ def parse_dict_with_units(cls, py_dict): my_dict_copy.pop("function") return super().parse_dict_with_units(my_dict_copy) - def operate(self, data): - return data - class Scale(abstractOperation): """ @@ -48,7 +45,8 @@ def operate(self, data, dep_var): """ Applies the operation for which the class is named for. - $f(\vec(x)) = scale*\vec(x)$ + .. math:: + f(\vec(x)) = scale*\vec(x) data: CSDM object dep_var: int. The index of the dependent variable to apply operation to @@ -62,7 +60,8 @@ class Gaussian(abstractOperation): Class for applying a Gaussian function to a dependent variable of simulation data - $f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2}$ + .. math:: + f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2} dimension: int. Data dimension to apply the function along sigma: float. Standard deviation of Gaussian function @@ -96,7 +95,8 @@ class Exponential(abstractOperation): Class for applying an exponential Lorentzian function to a dependent variable of simulation data - $f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)}$ + .. math:: + f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)} dimension: int. Data dimension to apply the function along Lambda: float. Width parameter diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 708e9bb1f..0eced7a61 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- -import numpy as np +import mrsimulator.post_simulation as ps from mrsimulator import Simulator +from mrsimulator.post_simulation import SignalProcessor try: from lmfit import Parameters @@ -128,19 +129,89 @@ def _traverse_dictionaries(dictionary, parent="spin_systems"): 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}")] - 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 + + 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(item.operations): + if isinstance(operation, ps.Gaussian): + identifier = prepend + f"opIndex_{i}_Gaussian" + arg = operation.sigma + temp_dict[f"{identifier}"] = arg + elif isinstance(operation, ps.Exponential): + identifier = prepend + f"opIndex_{i}_Exponential" + arg = operation.Lambda + temp_dict[f"{identifier}"] = arg + elif isinstance(operation, ps.Scale): + identifier = prepend + f"opIndex_{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_params(params, post_sim): + """ + Updates SignalProcessor operation arguments from an + LMFIT Parameters object + + params: LMFIT Parameters object + post_sim: SignalProcessor object + + """ + temp_dict = {} + arg_dict = {"Gaussian": "sigma", "Exponential": "Lambda", "Scale": "factor"} + for param in params: + # iterating through the parameter list looking for only DEP_VAR (ie post_sim params) + if "DEP_VAR" 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("opIndex") + 1] + val = params[param].value + # creating a dictionary of operations and arguments for each dependent variablle + if f"DepVar_{var}" not in temp_dict.keys(): + temp_dict[f"DepVar_{var}"] = {} + temp_dict[f"DepVar_{var}"][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[f"DepVar_{item.dependent_variable}"].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".operations[{split[0]}].{arg_dict[split[-1]]}" + assignment = f"={val}" + exec(dep_var_operation_list + operation_val_update + assignment) + + +def make_fitting_parameters(sim, post_sim=None, exclude_key=None): """ Parses through the fitting parameter list to create LMFIT parameters used for fitting. Args: sim: a Simulator object. + post_sim: a SignalProcessor object Returns: LMFIT Parameters object. @@ -155,69 +226,68 @@ def make_fitting_parameters(sim, exclude_key=None): if not isinstance(sim, Simulator): raise ValueError(f"Expecting a `Simulator` object, found {type(sim).__name__}.") + if not isinstance(post_sim, SignalProcessor) or post_sim is None: + raise ValueError( + f"Expecting a `SignalProcessor` object, found {type(post_sim).__name__}." + ) - 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")) - - 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")) - - for items in temp_list: - if "_eta" in items or "abundance" in items and last_abund not in items: - if "_eta" in items: + if isinstance(sim, Simulator): + params = Parameters() + temp_list = _traverse_dictionaries(_list_of_dictionaries(sim.spin_systems)) + + 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")) + + 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, + ) + elif last_abund 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, + expr=expression, ) - elif last_abund in items: - params.add( - name=items, - value=eval("sim." + _html_to_string(items)), - min=0, - max=100, - expr=expression, - ) - 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." + _html_to_string(items)) + if type(value) == list: + params.add(name=items, value=value[0]) + else: + params.add(name=items, value=value) + else: + params = Parameters() + temp_list = _traverse_dictionaries(_list_of_dictionaries(sim)) + + if isinstance(post_sim, SignalProcessor): + 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 min_function(params, sim, post_sim=None): """ The simulation routine to establish how the parameters will update the simulation. @@ -237,31 +307,44 @@ def min_function(params, sim, apodization_function=None): 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, 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: + if "DEP_VAR" not in items: nameString = "sim." + _html_to_string(items) executable = f"{nameString} = {values[items]}" exec(executable) - else: - nameString = "sim." + _html_to_string(items) - executable = f"{nameString} = [{values[items]}]" - exec(executable) + elif "DEP_VAR" in items and post_sim is not None: + _update_post_sim_from_params(params, post_sim) sim.run() - residual = np.asarray([]) + post_sim.data = sim.methods[0].simulation + post_sim.apply_operations() + # residual = np.asarray([]) + + if sim.config.decompose_spectrum == "spin_system": + for decomposed_datum in post_sim.data.dependent_variables: + datum = [sum(i) for i in zip(datum, decomposed_datum)] + else: + datum = post_sim.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 From 4e685a58a04e05a3f3e83516235fbf1a2cdafb8f Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Tue, 30 Jun 2020 11:48:26 -0400 Subject: [PATCH 07/26] Signal processing and associated documentation --- docs/api_py/py_api.rst | 1 + docs/api_py/signal_processing.rst | 15 ++ docs/load_spin_systems.rst | 4 +- .../plot_zpost_sim_signal_processing.py | 204 +++++++++++++++++ .../plot_1_mrsimFitExample_cuspidine.py | 43 ++-- .../Fitting/plot_2_mrsimFitExample_O17.py | 34 ++- src/mrsimulator/post_simulation.py | 211 ------------------ src/mrsimulator/signal_processing/__init__.py | 124 ++++++++++ src/mrsimulator/signal_processing/_base.py | 61 +++++ .../signal_processing/apodization.py | 123 ++++++++++ .../signal_processing}/test_apodization.py | 52 ++--- src/mrsimulator/simulator/__init__.py | 4 +- src/mrsimulator/spectral_fitting.py | 92 ++++---- tests/test_parameters.py | 42 ---- 14 files changed, 637 insertions(+), 373 deletions(-) create mode 100644 docs/api_py/signal_processing.rst create mode 100644 examples/1D_simulation/plot_zpost_sim_signal_processing.py delete mode 100644 src/mrsimulator/post_simulation.py create mode 100644 src/mrsimulator/signal_processing/__init__.py create mode 100644 src/mrsimulator/signal_processing/_base.py create mode 100644 src/mrsimulator/signal_processing/apodization.py rename {tests => src/mrsimulator/signal_processing}/test_apodization.py (59%) delete mode 100644 tests/test_parameters.py diff --git a/docs/api_py/py_api.rst b/docs/api_py/py_api.rst index e382ebe7a..8adf26a3e 100644 --- a/docs/api_py/py_api.rst +++ b/docs/api_py/py_api.rst @@ -15,6 +15,7 @@ Python-API References method methods fitting + signal_processing .. 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/load_spin_systems.rst b/docs/load_spin_systems.rst index 12e106c9c..c0f5c05d0 100644 --- a/docs/load_spin_systems.rst +++ b/docs/load_spin_systems.rst @@ -37,7 +37,7 @@ Here is an example. >>> from mrsimulator import Simulator >>> sim = Simulator() >>> filename = 'https://raw.githubusercontent.com/DeepanshS/mrsimulator-test/master/spin_systems_v0.3.json' - >>> sim.load_spin_systems(filename) + >>> sim.load_spin_systems(filename) # doctest:+SKIP Downloading '/DeepanshS/mrsimulator-test/master/spin_systems_v0.3.json' from 'raw.githubusercontent.com' to file 'spin_systems_v0.3.json'. [████████████████████████████████████] @@ -55,7 +55,7 @@ Here is an example. [████████████████████████████████████] >>> # The seven spin systems from the file are added to the sim object. - >>> len(sim.spin_systems) + >>> len(sim.spin_systems) # doctest:+SKIP 7 .. testsetup:: diff --git a/examples/1D_simulation/plot_zpost_sim_signal_processing.py b/examples/1D_simulation/plot_zpost_sim_signal_processing.py new file mode 100644 index 000000000..d8fc71c20 --- /dev/null +++ b/examples/1D_simulation/plot_zpost_sim_signal_processing.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Post Simulation Signal Processing. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. sectionauthor:: Maxwell C. Venetos +""" +# sphinx_gallery_thumbnail_number = 5 +# %% +# After running a simulation, we often want to process the resulting spectrum. +# For example, we may want to scale the intensities to match the experiment or +# apodize the signal to simulate line broadening. The following example will +# demonstrate the use of the `SignalProcessor` class to apply various operations +# to simulation data. +# global plot configuration +import matplotlib as mpl +import matplotlib.pyplot as plt + +font = {"weight": "light", "size": 9} +mpl.rc("font", **font) +mpl.rcParams["figure.figsize"] = [4.25, 3.0] + +# %% +# We will create a hypothetical two-site Si simulation to illustrate post-simulation +# signal processing. We will begin by processing the entire spectrum and follow up by +# decomposing the spectrum and processing each signal indepenently. + +from mrsimulator import SpinSystem +from mrsimulator import Simulator + +# %% +# **Step 1** Create the sites and spin system + +site1 = { + "isotope": "29Si", + "isotropic_chemical_shift": "-75.0 ppm", + "shielding_symmetric": {"zeta": "-60 ppm", "eta": 0.6}, +} + +site2 = { + "isotope": "29Si", + "isotropic_chemical_shift": "-80.0 ppm", + "shielding_symmetric": {"zeta": "-70 ppm", "eta": 0.5}, +} + +spin_system_1 = { + "name": "site A", + "description": "A test 29Si site", + "sites": [site1], # from the above code + "abundance": "100%", +} + + +spin_system_2 = { + "name": "site B", + "description": "A test 29Si site", + "sites": [site2], # from the above code + "abundance": "100%", +} + +system_object_1 = SpinSystem.parse_dict_with_units(spin_system_1) +system_object_2 = SpinSystem.parse_dict_with_units(spin_system_2) + +# %% +# **Step 2** Create the simulation and add the spin system objects. + +sim = Simulator() +sim.spin_systems += [system_object_1, system_object_2] + +# %% +# **Step 3** Create a Bloch decay spectrum method. + +from mrsimulator.methods import BlochDecaySpectrum + +method_dict = { + "channels": ["29Si"], + "magnetic_flux_density": "7.1 T", + "rotor_angle": "54.735 deg", + "rotor_frequency": "780 Hz", + "spectral_dimensions": [ + {"count": 2048, "spectral_width": "25 kHz", "reference_offset": "-5 kHz"} + ], +} +method_object = BlochDecaySpectrum.parse_dict_with_units(method_dict) +sim.methods += [method_object] + +# %% +# The above method is set up to record the :math:`^{29}\text{Si}` resonances at the +# magic angle, spinning at 780 Hz and 7.1 T external magnetic flux density. The +# resonances are recorded over 25 kHz spectral width using 2048 points. + +# %% +# **Step 4** Simulate the spectra. + +sim.run() + +# %% +# **Step 5** Plot the spectrum. + +ax = plt.subplot(projection="csdm") +ax.plot(sim.methods[0].simulation, color="black", linewidth=1) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + +# %% +# **Step 6** Create Post Simulation. + +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo + +# create list of processing operations +op_list1 = [ + sp.IFFT(dim_indx=0), + apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), +] + +post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list1) + +# %% +# The above signal processing procedure will apply the list of operations +# sequentially to the 0-index dependent variable to result in Gaussian +# broadening of the spectrum. The operation procedure will first perform +# a Fourier Transform to convert the data in frequency space to time space. +# Next, a Gaussian exponential with a broadening factor of 100 seconds is +# multipled to the data and then finally another Fourier transform is applied +# to the data to convert the data back to frequency space. + +# %% +# **Step 7** Apply the signal processing. + +processed_data = post_sim.apply_operations() + +# %% +# **Step 8** Plot the processed spectrum. + +ax = plt.subplot(projection="csdm") +ax.plot(processed_data, color="black", linewidth=1) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + +# %% +# The above code resulted in the same processing to be applied +# to both signals because in the simulation the signals were not +# seperated. 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() + +#%% +# **Step 5** Plot spectrum. + +x, y0, y1 = sim.methods[0].simulation.to_list() + +plt.plot(x, y0, x, y1) +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) +plt.tight_layout() +plt.show() + +#%% +# **Step 6** Create Post Simulation. + +op_list2 = [ + sp.IFFT(dim_indx=0), + apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), + apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=1), + sp.FFT(dim_indx=0), +] + +post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list2) + +# %% +# The above signal processing procedure will apply Gaussian +# broadening to dependent variable index 0 (-75 ppm Si site) +# and apply Lorentzian broadening to dependent variable +# index 1 (-80 ppm Si site). + +# %% +# **Step 7** Apply the singla processing. + +processed_data = post_sim.apply_operations() + +#%% +# **Step 8** Plot the processed spectrum + +x, y0, y1 = processed_data.to_list() + +plt.plot(x, y0, x, y1) +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) +plt.tight_layout() +plt.show() + +# %% diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index d21a9f502..560469f5e 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -109,29 +109,26 @@ # %% # **Step 6** Create a SignalProcessor -import mrsimulator.post_simulation as ps -from mrsimulator.post_simulation import SignalProcessor - -op_list = { - "dependent_variable": 0, - "operations": [ - ps.IFFT(dimension=0), - ps.Exponential(Lambda=200, dimension=0), - ps.FFT(dimension=0), - ps.Scale(factor=1), - ], -} +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo + +op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), + sp.Scale(factor=1), +] -post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) +post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) # %% # ** Step 7** Process and plot the spectrum. -post_sim.apply_operations() +processed_data = post_sim.apply_operations() ax = plt.subplot(projection="csdm") -ax.plot(post_sim.data, color="black", linewidth=1) +ax.plot(processed_data, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() @@ -164,10 +161,10 @@ # will showcase a fitting function provided in the *mrsimulator* library which automates the process. -def test_function(params, data, sim): +def test_function(params, sim, post_sim): values = params.valuesdict() - intensity = data.dependent_variables[0].components[0].real + 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] @@ -178,13 +175,12 @@ def test_function(params, data, sim): # here we run the simulation sim.run() - post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) - post_sim.operations[0].operations[3].factor = values["factor"] - post_sim.operations[0].operations[1].Lambda = values["Lambda"] + post_sim.operations[3].factor = values["factor"] + post_sim.operations[1].Lambda = values["Lambda"] # here we apodize the signal to simulate line broadening - post_sim.apply_operations() + processed_data = post_sim.apply_operations() - return intensity - post_sim.data.dependent_variables[0].components[0].real + return intensity - processed_data.dependent_variables[0].components[0].real #%% @@ -195,7 +191,7 @@ def test_function(params, data, sim): # %% **Step 8** Perform minimization. -minner = Minimizer(test_function, params, fcn_args=(synthetic_experiment, sim)) +minner = Minimizer(test_function, params, fcn_args=(sim, post_sim)) result = minner.minimize(method="powell") print(fit_report(result)) @@ -211,6 +207,7 @@ def test_function(params, data, sim): plt.plot(*residual.to_list(), alpha=0.5, label="Residual") plt.xlabel("Frequency / Hz") +plt.gca().invert_xaxis() plt.grid(which="major", axis="both", linestyle="--") plt.legend() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index a4d32b4e8..ab651c809 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -115,30 +115,27 @@ # %% # **Step 6** Create a SignalProcessor -import mrsimulator.post_simulation as ps -from mrsimulator.post_simulation import SignalProcessor +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo factor = oxygen_experiment.dependent_variables[0].components[0].max().real / 4 -op_list = { - "dependent_variable": 0, - "operations": [ - ps.IFFT(dimension=0), - ps.Exponential(Lambda=100, dimension=0), - ps.FFT(dimension=0), - ps.Scale(factor=factor), - ], -} +op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(Lambda=100, dim_indx=0, dep_var_indx=0), + sp.FFT(dim_indx=0), + sp.Scale(factor=factor), +] -post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[op_list]) +post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) # %% # ** Step 7** Process and plot the spectrum. -post_sim.apply_operations() +processed_data = post_sim.apply_operations() ax = plt.subplot(projection="csdm") -ax.plot(post_sim.data, color="black", linewidth=1) +ax.plot(processed_data, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() @@ -155,17 +152,17 @@ # **Step 8** Create a list of parameters to vary during fitting. -from mrsimulator.spectral_fitting import make_fitting_parameters +from mrsimulator.spectral_fitting import make_LMFIT_parameters -params = make_fitting_parameters(sim, post_sim) +params = make_LMFIT_parameters(sim, post_sim) params.pretty_print() # %% # **Step 9** Perform minimization. -from mrsimulator.spectral_fitting import min_function +from mrsimulator.spectral_fitting import LMFIT_min_function from lmfit import Minimizer, report_fit -minner = Minimizer(min_function, params, fcn_args=(sim, post_sim)) +minner = Minimizer(LMFIT_min_function, params, fcn_args=(sim, post_sim)) result = minner.minimize() report_fit(result) @@ -181,6 +178,7 @@ plt.plot(*residual.to_list(), alpha=0.5, label="Residual") plt.xlabel("$^{17}$O frequency / ppm") +plt.gca().invert_xaxis() plt.grid(which="major", axis="both", linestyle="--") plt.legend() diff --git a/src/mrsimulator/post_simulation.py b/src/mrsimulator/post_simulation.py deleted file mode 100644 index ad153b20d..000000000 --- a/src/mrsimulator/post_simulation.py +++ /dev/null @@ -1,211 +0,0 @@ -# -*- coding: utf-8 -*- -"""The Event class.""" -from typing import ClassVar -from typing import Dict -from typing import List -from typing import Union - -import csdmpy as cp -import numpy as np -from mrsimulator.util.parseable import Parseable -from pydantic import BaseModel - - -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.__class__.__name__ - return my_dict - - @classmethod - def parse_dict_with_units(cls, py_dict): - 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) - - -class Scale(abstractOperation): - """ - Class for applying a scaling - factor to a dependent variable of simulation data - """ - - factor: float = 1 - - def operate(self, data, dep_var): - """ - Applies the operation for which the class is named for. - - .. math:: - f(\vec(x)) = scale*\vec(x) - - data: CSDM object - dep_var: int. The index of the dependent variable to apply operation to - """ - data.dependent_variables[dep_var].components[0] *= self.factor - return data - - -class Gaussian(abstractOperation): - """ - Class for applying a Gaussian function - to a dependent variable of simulation data - - .. math:: - f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2} - - dimension: int. Data dimension to apply the function along - sigma: float. Standard deviation of Gaussian function - """ - - dimension: int = 0 - sigma: float = 0 - - property_unit_types: ClassVar = {"sigma": ["time", "frequency"]} - property_default_units: ClassVar = {"sigma": ["s", "Hz"]} - property_units: Dict = {"sigma": ["s", "Hz"]} - - def operate(self, data, dep_var): - """ - 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 - """ - # unit = self.property_units["sigma"] - # data_copy.dimensions[self.dimension].coordinates.to(unit).value - x = data.dimensions[self.dimension].coordinates.value - data.dependent_variables[dep_var].components[0] *= np.exp( - -((x * self.sigma * np.pi) ** 2) * 2 - ) - return data - - -class Exponential(abstractOperation): - """ - Class for applying an exponential - Lorentzian function to a dependent variable of simulation data - - .. math:: - f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)} - - dimension: int. Data dimension to apply the function along - Lambda: float. Width parameter - """ - - dimension: int = 0 - Lambda: float = 0 - - property_unit_types: ClassVar = {"Lambda": ["time", "frequency"]} - property_default_units: ClassVar = {"Lambda": ["s", "Hz"]} - property_units: Dict = {"Lambda": ["s", "Hz"]} - - def operate(self, data, dep_var): - """ - 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 - """ - # unit = self.property_units["Lambda"] - # x = data.dimensions[self.dimension].coordinates.to(unit).value - x = data.dimensions[self.dimension].coordinates.value - data.dependent_variables[dep_var].components[0] *= np.exp( - -self.Lambda * np.pi * np.abs(x) - ) - return data - - -class IFFT(abstractOperation): - """ - Class for applying an inverse - Fourier transform to a dependent variable of simulation data - - dimension: int. Data dimension to apply the function along - - """ - - dimension: int = 0 - - def operate(self, data, dep_var): - """ - 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 data.fft(axis=self.dimension) - - -class FFT(abstractOperation): - """ - Class for applying a - Fourier transform to a dependent variable of simulation data - - dimension: int. Data dimension to apply the function along - - """ - - dimension: int = 0 - - def operate(self, data, dep_var): - """ - 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 data.fft(axis=self.dimension) - - -class OperationList(BaseModel): - """ - Container class to hold a list of operations to be applied to a given - dependent variable in the data - - dependent_variable: The dependent variable of the data to apply operations to - operations: List of abstractOperation based operations to sequentially apply to data. - """ - - dependent_variable: int = 0 - operations: List[abstractOperation] - - -class SignalProcessor(BaseModel): - """ - Signal processing class to apply lists of various operations to individual dependent variables - of the data. - - data: CSDM object. From simulation - operations: List of operation lists - """ - - data: Union[cp.CSDM, np.ndarray] = None - operations: List[OperationList] - - class Config: - validate_assignment = True - arbitrary_types_allowed = True - - def apply_operations(self, **kwargs): - op_list = self.operations - - copy_data = self.data.copy() - - for item in op_list: - dep_var = item.dependent_variable - for filters in item.operations: - copy_data = filters.operate(copy_data, dep_var=dep_var) - - self.data = copy_data - - return copy_data diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py new file mode 100644 index 000000000..36d327a08 --- /dev/null +++ b/src/mrsimulator/signal_processing/__init__.py @@ -0,0 +1,124 @@ +# -*- coding: utf-8 -*- +"""The Event class.""" +from typing import List +from typing import Union + +import csdmpy as cp +import numpy as np +from pydantic import BaseModel + +from ._base import abstractOperation + +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" + + +class SignalProcessor(BaseModel): + """ + Signal processing class to apply lists of various operations to individual dependent variables + of the data. + + Attributes + ---------- + + data: CSDM object. + From simulation + + operations: List + List of operation lists + + Examples + -------- + + >>> post_sim = SignalProcessor(data = csdm_object, operations = [op1, op2]) # doctest: +SKIP + """ + + data: Union[cp.CSDM, np.ndarray] = None + operations: List[abstractOperation] + + class Config: + validate_assignment = True + arbitrary_types_allowed = True + + 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, **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. + """ + copy_data = self.data.copy() + for filters in self.operations: + copy_data = filters.operate(copy_data) + + return copy_data + + +class Scale(abstractOperation): + """ + Class for applying a scaling + factor to a dependent variable of simulation data + """ + + 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): + """ + Class for applying an inverse + Fourier transform to a dependent variable of simulation data + + Args: + dim_indx: int. Data dimension to apply the function along + + """ + + 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): + 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..5e0753afe --- /dev/null +++ b/src/mrsimulator/signal_processing/_base.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- +"""The abstractOperation class.""" +from mrsimulator.util.parseable import Parseable + + +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): + 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) + + +# class FFT(abstractOperation): +# """ +# Class for applying a +# Fourier transform to a dependent variable of simulation data + +# dimension: int. Data dimension to apply the function along + +# """ + +# dimension: int = 0 + +# 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 data.fft(axis=self.dimension) + + +# class OperationList(BaseModel): +# """ +# Container class to hold a list of operations to be applied to a given +# dependent variable in the data + +# dependent_variable: The dependent variable of the data to apply operations to +# operations: List of abstractOperation based operations to sequentially apply to data. +# """ + +# dependent_variable: int = 0 +# operations: List[abstractOperation] diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py new file mode 100644 index 000000000..70fdac9dc --- /dev/null +++ b/src/mrsimulator/signal_processing/apodization.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +"""The Event class.""" +from typing import ClassVar +from typing import Dict + +import numpy as np + +from ._base import abstractOperation + + +class AbstractApodization(abstractOperation): + @property + def function(self): + return "apodization" + + @property + def type(self): + return self.__class__.__name__ + + +class Gaussian(AbstractApodization): + r""" + Class for applying a Gaussian function + to a dependent variable of simulation data + + .. math:: + f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2} + + Args: + dim_indx: int. Data dimension to apply the function along. + sigma: float. Standard deviation of Gaussian function + dep_var_indx: int. Data dependent variable index to apply the function to. + If type None, will be applied to every dependent variable. + + """ + + dim_indx: int = 0 + sigma: float = 0 + dep_var_indx: int = None # if none apply to all + + property_unit_types: ClassVar = {"sigma": ["time", "frequency"]} + property_default_units: ClassVar = {"sigma": ["s", "Hz"]} + property_units: Dict = {"sigma": "Hz"} + + 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 + """ + # unit = self.property_units["sigma"] + # data_copy.dimensions[self.dimension].coordinates.to(unit).value + x = data.dimensions[self.dim_indx].coordinates + x_value, x_unit = x.value, x.unit + if x_unit.physical_type == "frequency": + self.property_units = dict(sigma="s") + else: + self.property_units = dict(sigma="Hz") + + if self.dep_var_indx is None: + for i in range(len(data.dependent_variables)): + data.dependent_variables[i].components[0] *= np.exp( + -((x_value * self.sigma * np.pi) ** 2) * 2 + ) + + else: + data.dependent_variables[self.dep_var_indx].components[0] *= np.exp( + -((x_value * self.sigma * np.pi) ** 2) * 2 + ) + return data + + +class Exponential(AbstractApodization): + r""" + Class for applying an exponential + Lorentzian function to a dependent variable of simulation data + + .. math:: + f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)} + + Args: + dim_indx: int. Data dimension to apply the function along. + Lambda: float. Width parameter + dep_var_indx: int. Data dependent variable index to apply the function to. + If type None, will be applied to every dependent variable. + + """ + + dim_indx: int = 0 + Lambda: float = 0 + dep_var_indx: int = None # if none apply to all + + property_unit_types: ClassVar = {"Lambda": ["time", "frequency"]} + property_default_units: ClassVar = {"Lambda": ["s", "Hz"]} + property_units: Dict = {"Lambda": "Hz"} + + 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 + """ + # unit = self.property_units["Lambda"] + # x = data.dimensions[self.dimension].coordinates.to(unit).value + x = data.dimensions[self.dim_indx].coordinates + x_value, x_unit = x.value, x.unit + if x_unit.physical_type == "frequency": + self.property_units = dict(Lambda="s") + else: + self.property_units = dict(Lambda="Hz") + + if self.dep_var_indx is None: + for i in range(len(data.dependent_variables)): + data.dependent_variables[i].components[0] *= np.exp( + -self.Lambda * np.pi * np.abs(x_value) + ) + else: + data.dependent_variables[self.dep_var_indx].components[0] *= np.exp( + -self.Lambda * np.pi * np.abs(x_value) + ) + return data diff --git a/tests/test_apodization.py b/src/mrsimulator/signal_processing/test_apodization.py similarity index 59% rename from tests/test_apodization.py rename to src/mrsimulator/signal_processing/test_apodization.py index a49f88de6..c35b1a85c 100644 --- a/tests/test_apodization.py +++ b/src/mrsimulator/signal_processing/test_apodization.py @@ -1,11 +1,11 @@ # -*- coding: utf-8 -*- """Apodization test""" -import mrsimulator.post_simulation as ps +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 -from mrsimulator.post_simulation import SignalProcessor sim = Simulator() the_site = {"isotope": "1H", "isotropic_chemical_shift": "0 ppm"} @@ -23,25 +23,19 @@ ], ) -PS_0 = {"dependent_variable": 0, "operations": [ps.Scale(factor=10)]} +PS_0 = [sp.Scale(factor=10)] -PS_1 = { - "dependent_variable": 0, - "operations": [ - ps.IFFT(dimension=0), - ps.Exponential(Lambda=200, dimension=0), - ps.FFT(dimension=0), - ], -} - -PS_2 = { - "dependent_variable": 0, - "operations": [ - ps.IFFT(dimension=0), - ps.Gaussian(sigma=20, dimension=0), - ps.FFT(dimension=0), - ], -} +PS_1 = [ + sp.IFFT(dim_indx=0), + apo.Exponential(Lambda=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), + sp.FFT(dim_indx=0), +] sim.methods += [method_1] sim.run() @@ -51,18 +45,18 @@ def test_scale(): - post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_0]) - post_sim.apply_operations() + post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_0) + data = post_sim.apply_operations() x0, y0 = sim.methods[0].simulation.to_list() - x, y = post_sim.data.to_list() + x, y = data.to_list() assert y.max() / y0.max() == 10, "Scaling failed" def test_Lorentzian(): - post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_1]) - post_sim.apply_operations() - x, y = post_sim.data.to_list() + post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_1) + data = post_sim.apply_operations() + x, y = data.to_list() sigma = 200 test = (sigma / 2) / (np.pi * (freqHz ** 2 + (sigma / 2) ** 2)) @@ -73,9 +67,9 @@ def test_Lorentzian(): def test_Gaussian(): - post_sim = SignalProcessor(data=sim.methods[0].simulation, operations=[PS_2]) - post_sim.apply_operations() - x, y = post_sim.data.to_list() + post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_2) + data = post_sim.apply_operations() + x, y = data.to_list() sigma = 20 test = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((freqHz / sigma) ** 2) / 2) diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index dcbd20410..07a8d6ec4 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -324,8 +324,8 @@ def load_spin_systems(self, filename: str): >>> sim.load_spin_systems(filename) # doctest:+SKIP """ contents = import_json(filename) - json_data = contents["spin_systems"] - self.spin_systems = [SpinSystem.parse_dict_with_units(obj) for obj in json_data] + # json_data = contents["spin_systems"] + self.spin_systems = [SpinSystem.parse_dict_with_units(obj) for obj in contents] def export_spin_systems(self, filename: str): """ diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 0eced7a61..c8c51deb1 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- -import mrsimulator.post_simulation as ps +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo from mrsimulator import Simulator -from mrsimulator.post_simulation import SignalProcessor try: from lmfit import Parameters @@ -142,21 +142,21 @@ def _post_sim_LMFIT_params(post_sim): returns: Parameters object """ temp_dict = {} - for item in post_sim.operations: - prepend = f"DEP_VAR_{item.dependent_variable}_" - for i, operation in enumerate(item.operations): - if isinstance(operation, ps.Gaussian): - identifier = prepend + f"opIndex_{i}_Gaussian" - arg = operation.sigma - temp_dict[f"{identifier}"] = arg - elif isinstance(operation, ps.Exponential): - identifier = prepend + f"opIndex_{i}_Exponential" - arg = operation.Lambda - temp_dict[f"{identifier}"] = arg - elif isinstance(operation, ps.Scale): - identifier = prepend + f"opIndex_{i}_Scale" - arg = operation.factor - temp_dict[f"{identifier}"] = arg + # 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.Lambda + 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(): @@ -165,7 +165,7 @@ def _post_sim_LMFIT_params(post_sim): return params -def _update_post_sim_from_params(params, post_sim): +def _update_post_sim_from_LMFIT_params(params, post_sim): """ Updates SignalProcessor operation arguments from an LMFIT Parameters object @@ -178,33 +178,33 @@ def _update_post_sim_from_params(params, post_sim): arg_dict = {"Gaussian": "sigma", "Exponential": "Lambda", "Scale": "factor"} for param in params: # iterating through the parameter list looking for only DEP_VAR (ie post_sim params) - if "DEP_VAR" in param: + 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("opIndex") + 1] + # 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 variablle - if f"DepVar_{var}" not in temp_dict.keys(): - temp_dict[f"DepVar_{var}"] = {} - temp_dict[f"DepVar_{var}"][f"{opIndex}_{split_name[-1]}"] = val + # 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[f"DepVar_{item.dependent_variable}"].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".operations[{split[0]}].{arg_dict[split[-1]]}" - assignment = f"={val}" - exec(dep_var_operation_list + operation_val_update + assignment) - - -def make_fitting_parameters(sim, post_sim=None, exclude_key=None): + # 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 through the fitting parameter list to create LMFIT parameters used for fitting. @@ -226,7 +226,7 @@ def make_fitting_parameters(sim, post_sim=None, exclude_key=None): if not isinstance(sim, Simulator): raise ValueError(f"Expecting a `Simulator` object, found {type(sim).__name__}.") - if not isinstance(post_sim, SignalProcessor) or post_sim is None: + if not isinstance(post_sim, sp.SignalProcessor) or post_sim is None: raise ValueError( f"Expecting a `SignalProcessor` object, found {type(post_sim).__name__}." ) @@ -278,7 +278,7 @@ def make_fitting_parameters(sim, post_sim=None, exclude_key=None): params = Parameters() temp_list = _traverse_dictionaries(_list_of_dictionaries(sim)) - if isinstance(post_sim, SignalProcessor): + if isinstance(post_sim, sp.SignalProcessor): temp_params = _post_sim_LMFIT_params(post_sim) for item in temp_params: params.add(name=item, value=temp_params[item].value) @@ -287,7 +287,7 @@ def make_fitting_parameters(sim, post_sim=None, exclude_key=None): return params -def min_function(params, sim, post_sim=None): +def LMFIT_min_function(params, sim, post_sim=None): """ The simulation routine to establish how the parameters will update the simulation. @@ -307,7 +307,7 @@ def min_function(params, sim, post_sim=None): raise ValueError( f"Expecting a `Parameters` object, found {type(params).__name__}." ) - if not isinstance(post_sim, SignalProcessor) or post_sim is None: + if not isinstance(post_sim, sp.SignalProcessor) or post_sim is None: raise ValueError( f"Expecting a `SignalProcessor` object, found {type(post_sim).__name__}." ) @@ -317,23 +317,23 @@ def min_function(params, sim, post_sim=None): values = params.valuesdict() for items in values: - if "DEP_VAR" not in items: + if "operation_" not in items: nameString = "sim." + _html_to_string(items) executable = f"{nameString} = {values[items]}" exec(executable) - elif "DEP_VAR" in items and post_sim is not None: - _update_post_sim_from_params(params, post_sim) + elif "operation_" in items and post_sim is not None: + _update_post_sim_from_LMFIT_params(params, post_sim) sim.run() post_sim.data = sim.methods[0].simulation - post_sim.apply_operations() + processed_data = post_sim.apply_operations() # residual = np.asarray([]) if sim.config.decompose_spectrum == "spin_system": - for decomposed_datum in post_sim.data.dependent_variables: + for decomposed_datum in processed_data.dependent_variables: datum = [sum(i) for i in zip(datum, decomposed_datum)] else: - datum = post_sim.data.dependent_variables[0].components[0] + datum = processed_data.dependent_variables[0].components[0] return ( sim.methods[0].experiment.dependent_variables[0].components[0].real - datum.real 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" From b9785cefb3ffb98cc95ad37a53037ccc8327cac1 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 5 Jul 2020 22:12:11 -0400 Subject: [PATCH 08/26] fixed pytest failures. --- .flake8 | 2 +- docs/load_spin_systems.rst | 8 ++++---- src/mrsimulator/method/__init__.py | 4 ++-- src/mrsimulator/spectral_fitting.py | 6 ++++-- 4 files changed, 11 insertions(+), 9 deletions(-) diff --git a/.flake8 b/.flake8 index e598e9d8e..6c322c13f 100644 --- a/.flake8 +++ b/.flake8 @@ -3,5 +3,5 @@ ignore = E402 #C901 E265 E402 E501 exclude = .eggs, *.egg,build, src/mrsimulator/__init__.py filename = *.pyx, *py max-line-length = 88 -max-complexity = 10 +max-complexity = 18 select = C,E,F,W,N8 diff --git a/docs/load_spin_systems.rst b/docs/load_spin_systems.rst index c0f5c05d0..1834be17f 100644 --- a/docs/load_spin_systems.rst +++ b/docs/load_spin_systems.rst @@ -37,10 +37,10 @@ Here is an example. >>> from mrsimulator import Simulator >>> sim = Simulator() >>> filename = 'https://raw.githubusercontent.com/DeepanshS/mrsimulator-test/master/spin_systems_v0.3.json' - >>> sim.load_spin_systems(filename) # doctest:+SKIP + >>> sim.load_spin_systems(filename) Downloading '/DeepanshS/mrsimulator-test/master/spin_systems_v0.3.json' from 'raw.githubusercontent.com' to file 'spin_systems_v0.3.json'. - [████████████████████████████████████] + [███████████████████████] .. doctest:: @@ -52,10 +52,10 @@ Here is an example. >>> sim.load_spin_systems(filename) # doctest:+SKIP Downloading '/DeepanshS/mrsimulator-test/master/spin_systems_v0.3.json' from 'raw.githubusercontent.com' to file 'spin_systems_v0.3.json'. - [████████████████████████████████████] + [███████████████████████] >>> # The seven spin systems from the file are added to the sim object. - >>> len(sim.spin_systems) # doctest:+SKIP + >>> len(sim.spin_systems) 7 .. testsetup:: diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 87804d268..f394d1d8d 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -41,8 +41,8 @@ class Method(Parseable): Example ------- - >>> bloch = Method() - >>> bloch.channels = ['1H'] + >>> bloch = Method() + >>> bloch.channels = ['1H'] spectral_dimensions: list of :ref:`spectral_dim_api` or dict objects (optional). The number of spectral dimensions depends on the given method. For example, a diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index c8c51deb1..37eac9e32 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -187,7 +187,7 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): # 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 variablle + # 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 @@ -330,8 +330,10 @@ def LMFIT_min_function(params, sim, post_sim=None): # residual = np.asarray([]) if sim.config.decompose_spectrum == "spin_system": + datum = 0 for decomposed_datum in processed_data.dependent_variables: - datum = [sum(i) for i in zip(datum, decomposed_datum)] + datum += decomposed_datum.components[0] + # datum = [sum(i) for i in zip(datum, decomposed_datum)] else: datum = processed_data.dependent_variables[0].components[0] From d98721b6573ca812dc99f118aea21d9d95c5a1d0 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 5 Jul 2020 22:18:25 -0400 Subject: [PATCH 09/26] simplyfied _str_to_html and _html_to_str. --- src/mrsimulator/spectral_fitting.py | 69 ++++++++++++----------------- 1 file changed, 29 insertions(+), 40 deletions(-) diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 37eac9e32..ebf154a4d 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -15,6 +15,30 @@ __email__ = "maxvenetos@gmail.com" +ENCRYPTION_PAIRS = [ + ["spin_systems[", "ISO_"], + ["].sites[", "_SITES_"], + ["].isotropic_chemical_shift", "_isotropic_chemical_shift"], + ["].shielding_symmetric.", "_shielding_symmetric_"], + ["].quadrupolar.", "_quadrupolar_"], + ["].abundance", "_abundance"], + ["methods[", "METHODS_"], # why does methods needs 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): """ LMFIT Parameters class does not allow for names to include special characters. @@ -28,20 +52,8 @@ def _str_to_html(my_string): 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 ENCRYPTION_PAIRS: + my_string = my_string.replace(*item) return my_string @@ -57,21 +69,8 @@ def _html_to_string(my_string): 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 ENCRYPTION_PAIRS: + my_string = my_string.replace(*item[::-1]) return my_string @@ -90,16 +89,6 @@ def _list_of_dictionaries(my_list): 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 @@ -118,7 +107,7 @@ def _traverse_dictionaries(dictionary, parent="spin_systems"): 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}") From 89142ddb045577bd2d230a1b0b79b6998f777686 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 5 Jul 2020 22:53:46 -0400 Subject: [PATCH 10/26] simplified make_LMFIT_parameters function. --- .flake8 | 2 +- src/mrsimulator/spectral_fitting.py | 105 +++++++++++++--------------- 2 files changed, 49 insertions(+), 58 deletions(-) diff --git a/.flake8 b/.flake8 index 6c322c13f..d7840d4ad 100644 --- a/.flake8 +++ b/.flake8 @@ -3,5 +3,5 @@ ignore = E402 #C901 E265 E402 E501 exclude = .eggs, *.egg,build, src/mrsimulator/__init__.py filename = *.pyx, *py max-line-length = 88 -max-complexity = 18 +max-complexity = 12 select = C,E,F,W,N8 diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index ebf154a4d..13ff24790 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -14,9 +14,9 @@ __author__ = "Maxwell C Venetos" __email__ = "maxvenetos@gmail.com" - +START = "ISO_" ENCRYPTION_PAIRS = [ - ["spin_systems[", "ISO_"], + ["spin_systems[", START], ["].sites[", "_SITES_"], ["].isotropic_chemical_shift", "_isotropic_chemical_shift"], ["].shielding_symmetric.", "_shielding_symmetric_"], @@ -209,69 +209,60 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): 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) if not isinstance(sim, Simulator): raise ValueError(f"Expecting a `Simulator` object, found {type(sim).__name__}.") - if not isinstance(post_sim, sp.SignalProcessor) or post_sim is None: + + params = Parameters() + temp_list = _traverse_dictionaries(_list_of_dictionaries(sim.spin_systems)) + + # get total abundance scaling factor + length = len(sim.spin_systems) + 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 = "100" + "-".join([f"{START}{i}_abundance" for i in range(length - 1)]) + + for items in temp_list: + if "_eta" in items: + params.add( + name=items, value=eval("sim." + _html_to_string(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)), + min=0, + max=100, + expr=expression, + ) + elif "abundance" in items: + params.add( + name=items, + value=eval("sim." + _html_to_string(items)) * abundance_scale, + min=0, + max=100, + ) + else: + value = eval("sim." + _html_to_string(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__}." ) - - if isinstance(sim, Simulator): - params = Parameters() - temp_list = _traverse_dictionaries(_list_of_dictionaries(sim.spin_systems)) - - 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")) - - 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, - ) - elif last_abund in items: - params.add( - name=items, - value=eval("sim." + _html_to_string(items)), - min=0, - max=100, - expr=expression, - ) - 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) - else: - params = Parameters() - temp_list = _traverse_dictionaries(_list_of_dictionaries(sim)) - - if isinstance(post_sim, sp.SignalProcessor): - 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) + 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 From 0f25ba49c6bcc6e3a632f9022e402592aadb64ad Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 5 Jul 2020 23:45:21 -0400 Subject: [PATCH 11/26] Fixed python api docs for fitting utility. --- docs/api_py/fitting.rst | 4 ++-- docs/api_py/simulator.rst | 1 - src/mrsimulator/signal_processing/_base.py | 5 +++++ 3 files changed, 7 insertions(+), 3 deletions(-) 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/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/src/mrsimulator/signal_processing/_base.py b/src/mrsimulator/signal_processing/_base.py index 5e0753afe..cf658c205 100644 --- a/src/mrsimulator/signal_processing/_base.py +++ b/src/mrsimulator/signal_processing/_base.py @@ -21,6 +21,11 @@ def to_dict_with_units(self): @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") From 5c5158b249d3128db16f88bf528aeca673cad225 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 6 Jul 2020 00:08:08 -0400 Subject: [PATCH 12/26] Rename Zeeman_state.py to zeeman_state.py --- src/mrsimulator/spin_system/{Zeeman_state.py => zeeman_state.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/mrsimulator/spin_system/{Zeeman_state.py => zeeman_state.py} (100%) diff --git a/src/mrsimulator/spin_system/Zeeman_state.py b/src/mrsimulator/spin_system/zeeman_state.py similarity index 100% rename from src/mrsimulator/spin_system/Zeeman_state.py rename to src/mrsimulator/spin_system/zeeman_state.py From 25e19c96699416253261b3a9ce132da6c2f4bdf5 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 6 Jul 2020 14:51:20 -0400 Subject: [PATCH 13/26] add apodization class tests simplified apodization class. --- src/mrsimulator/signal_processing/__init__.py | 4 +- .../signal_processing/apodization.py | 120 +++++++----- .../signal_processing/test_apodization.py | 79 -------- .../tests/test_apodization.py | 178 ++++++++++++++++++ 4 files changed, 256 insertions(+), 125 deletions(-) delete mode 100644 src/mrsimulator/signal_processing/test_apodization.py create mode 100644 src/mrsimulator/signal_processing/tests/test_apodization.py diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py index 36d327a08..f6179caa5 100644 --- a/src/mrsimulator/signal_processing/__init__.py +++ b/src/mrsimulator/signal_processing/__init__.py @@ -1,10 +1,8 @@ # -*- coding: utf-8 -*- """The Event class.""" from typing import List -from typing import Union import csdmpy as cp -import numpy as np from pydantic import BaseModel from ._base import abstractOperation @@ -33,7 +31,7 @@ class SignalProcessor(BaseModel): >>> post_sim = SignalProcessor(data = csdm_object, operations = [op1, op2]) # doctest: +SKIP """ - data: Union[cp.CSDM, np.ndarray] = None + data: cp.CSDM = None operations: List[abstractOperation] class Config: diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index 70fdac9dc..2e9ec066e 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -1,7 +1,10 @@ # -*- coding: utf-8 -*- """The Event class.""" +from copy import deepcopy +from sys import modules from typing import ClassVar from typing import Dict +from typing import Union import numpy as np @@ -9,6 +12,17 @@ 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): + copy_dict = deepcopy(py_dict) + copy_dict.pop("type") + + obj = super().parse_dict_with_units(copy_dict) + return getattr(modules[__name__], py_dict["type"])(**obj.dict()) + @property def function(self): return "apodization" @@ -17,6 +31,53 @@ def function(self): def type(self): 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 = dict(prop="s") + return + self.property_units = dict(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""" @@ -34,14 +95,16 @@ class Gaussian(AbstractApodization): """ - dim_indx: int = 0 sigma: float = 0 - dep_var_indx: int = None # if none apply to all 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(-((x * arg * np.pi) ** 2) * 2) + def operate(self, data): """ Applies the operation for which the class is named for. @@ -49,26 +112,10 @@ def operate(self, data): data: CSDM object dep_var: int. The index of the dependent variable to apply operation to """ - # unit = self.property_units["sigma"] - # data_copy.dimensions[self.dimension].coordinates.to(unit).value - x = data.dimensions[self.dim_indx].coordinates - x_value, x_unit = x.value, x.unit - if x_unit.physical_type == "frequency": - self.property_units = dict(sigma="s") - else: - self.property_units = dict(sigma="Hz") - - if self.dep_var_indx is None: - for i in range(len(data.dependent_variables)): - data.dependent_variables[i].components[0] *= np.exp( - -((x_value * self.sigma * np.pi) ** 2) * 2 - ) - - else: - data.dependent_variables[self.dep_var_indx].components[0] *= np.exp( - -((x_value * self.sigma * np.pi) ** 2) * 2 - ) - return data + + return super().operate( + data, fn=self.fn, prop_name="sigma", prop_value=self.sigma + ) class Exponential(AbstractApodization): @@ -87,14 +134,16 @@ class Exponential(AbstractApodization): """ - dim_indx: int = 0 Lambda: float = 0 - dep_var_indx: int = None # if none apply to all property_unit_types: ClassVar = {"Lambda": ["time", "frequency"]} property_default_units: ClassVar = {"Lambda": ["s", "Hz"]} property_units: Dict = {"Lambda": "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. @@ -102,22 +151,7 @@ def operate(self, data): data: CSDM object dep_var: int. The index of the dependent variable to apply operation to """ - # unit = self.property_units["Lambda"] - # x = data.dimensions[self.dimension].coordinates.to(unit).value - x = data.dimensions[self.dim_indx].coordinates - x_value, x_unit = x.value, x.unit - if x_unit.physical_type == "frequency": - self.property_units = dict(Lambda="s") - else: - self.property_units = dict(Lambda="Hz") - - if self.dep_var_indx is None: - for i in range(len(data.dependent_variables)): - data.dependent_variables[i].components[0] *= np.exp( - -self.Lambda * np.pi * np.abs(x_value) - ) - else: - data.dependent_variables[self.dep_var_indx].components[0] *= np.exp( - -self.Lambda * np.pi * np.abs(x_value) - ) - return data + + return super().operate( + data, fn=self.fn, prop_name="Lambda", prop_value=self.Lambda + ) diff --git a/src/mrsimulator/signal_processing/test_apodization.py b/src/mrsimulator/signal_processing/test_apodization.py deleted file mode 100644 index c35b1a85c..000000000 --- a/src/mrsimulator/signal_processing/test_apodization.py +++ /dev/null @@ -1,79 +0,0 @@ -# -*- 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] - -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(Lambda=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), - 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(data=sim.methods[0].simulation, operations=PS_0) - data = post_sim.apply_operations() - x0, y0 = sim.methods[0].simulation.to_list() - x, y = data.to_list() - - assert y.max() / y0.max() == 10, "Scaling failed" - - -def test_Lorentzian(): - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_1) - data = post_sim.apply_operations() - x, y = data.to_list() - - sigma = 200 - test = (sigma / 2) / (np.pi * (freqHz ** 2 + (sigma / 2) ** 2)) - - assert np.allclose( - test / test.max(), y / y.max(), atol=1e-04 - ), "Lorentzian apodization amplitude failed" - - -def test_Gaussian(): - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_2) - data = post_sim.apply_operations() - x, y = data.to_list() - - sigma = 20 - test = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-((freqHz / sigma) ** 2) / 2) - - assert np.allclose( - test / test.max(), y / y.max(), atol=1e-04 - ), "Gaussian apodization amplitude failed" 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..78bf9ac99 --- /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(Lambda=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(data=sim.methods[0].simulation, operations=PS_0) + data = post_sim.apply_operations() + _, 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(data=sim.methods[0].simulation, operations=PS_1) + data = post_sim.apply_operations() + _, 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(data=sim.methods[0].simulation, operations=PS_2) + data = post_sim.apply_operations() + _, 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(data=sim.methods[0].simulation, operations=PS_3) + data = post_sim.apply_operations() + _, 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(Lambda=200, dim_indx=0, dep_var_indx=0) + + assert a.Lambda == 200 + assert a.property_units == {"Lambda": "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", + "Lambda": "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 From 396c77a9c45954f68d29533c594cf3148bfb6fa5 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 6 Jul 2020 19:33:34 -0400 Subject: [PATCH 14/26] add SignalProcessor class tests. minor code style formatting. --- .../plot_zpost_sim_signal_processing.py | 6 +-- .../plot_1_mrsimFitExample_cuspidine.py | 12 ++--- .../Fitting/plot_2_mrsimFitExample_O17.py | 10 ++-- src/mrsimulator/signal_processing/__init__.py | 51 ++++++++++++------- src/mrsimulator/signal_processing/_base.py | 46 +++-------------- .../signal_processing/apodization.py | 26 ++++------ .../tests/test_signal_processing.py | 35 +++++++++++++ src/mrsimulator/simulator/__init__.py | 2 +- src/mrsimulator/spectral_fitting.py | 2 +- src/mrsimulator/spin_system/__init__.py | 3 +- src/mrsimulator/spin_system/tensors.py | 2 +- 11 files changed, 106 insertions(+), 89 deletions(-) create mode 100644 src/mrsimulator/signal_processing/tests/test_signal_processing.py diff --git a/examples/1D_simulation/plot_zpost_sim_signal_processing.py b/examples/1D_simulation/plot_zpost_sim_signal_processing.py index d8fc71c20..04fd5b569 100644 --- a/examples/1D_simulation/plot_zpost_sim_signal_processing.py +++ b/examples/1D_simulation/plot_zpost_sim_signal_processing.py @@ -154,7 +154,7 @@ sim.config.decompose_spectrum = "spin_system" sim.run() -#%% +# %% # **Step 5** Plot spectrum. x, y0, y1 = sim.methods[0].simulation.to_list() @@ -166,7 +166,7 @@ plt.tight_layout() plt.show() -#%% +# %% # **Step 6** Create Post Simulation. op_list2 = [ @@ -189,7 +189,7 @@ processed_data = post_sim.apply_operations() -#%% +# %% # **Step 8** Plot the processed spectrum x, y0, y1 = processed_data.to_list() diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 560469f5e..d6608daac 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -133,7 +133,7 @@ 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 @@ -141,7 +141,7 @@ # examples on fitting. Here, however, we will construct the parameter list explicitly to demonstrate how the parameters # are created. -#%% +# %% from lmfit import Minimizer, Parameters, fit_report @@ -155,7 +155,7 @@ params.add(name="factor", value=1) params -#%% +# %% # 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. @@ -183,7 +183,7 @@ def test_function(params, sim, post_sim): return intensity - processed_data.dependent_variables[0].components[0].real -#%% +# %% # 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 magic angle spinning fitting is we must use a discrete minimization method @@ -195,7 +195,7 @@ def test_function(params, sim, post_sim): result = minner.minimize(method="powell") print(fit_report(result)) -#%% +# %% # **Step 9** Plot the fitted spectrum. @@ -214,7 +214,7 @@ def test_function(params, sim, post_sim): 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 diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index ab651c809..489a1901a 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -6,7 +6,7 @@ .. 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 @@ -26,7 +26,7 @@ mpl.rc("font", **font) mpl.rcParams["figure.figsize"] = [4.25, 3.0] -#%% +# %% # Next we will import `csdmpy `_ and loading the data file. import csdmpy as cp @@ -42,7 +42,7 @@ 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. @@ -140,7 +140,7 @@ 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*. # @@ -185,7 +185,7 @@ 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/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py index f6179caa5..1bebf0887 100644 --- a/src/mrsimulator/signal_processing/__init__.py +++ b/src/mrsimulator/signal_processing/__init__.py @@ -1,11 +1,12 @@ # -*- 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 +from ._base import AbstractOperation __author__ = "Maxwell C. Venetos" __email__ = "maxvenetos@gmail.com" @@ -13,8 +14,8 @@ class SignalProcessor(BaseModel): """ - Signal processing class to apply lists of various operations to individual dependent variables - of the data. + Signal processing class to apply lists of various operations to individual + dependent variables of the data. Attributes ---------- @@ -32,12 +33,33 @@ class SignalProcessor(BaseModel): """ data: cp.CSDM = None - operations: List[abstractOperation] + 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 @@ -69,17 +91,15 @@ def apply_operations(self, **kwargs): return copy_data -class Scale(abstractOperation): +class Scale(AbstractOperation): """ - Class for applying a scaling - factor to a dependent variable of simulation data + Class for applying a scaling factor to a dependent variable of simulation data. """ factor: float = 1 def operate(self, data): - r""" - Applies the operation for which the class is named for. + r"""Applies the operation for which the class is named for. .. math:: f(\vec(x)) = scale*\vec(x) @@ -91,25 +111,22 @@ def operate(self, data): return data -class IFFT(abstractOperation): +class IFFT(AbstractOperation): """ - Class for applying an inverse - Fourier transform to a dependent variable of simulation data + Class for applying an inverse Fourier transform to a dependent variable of + simulation data. Args: dim_indx: int. Data dimension to apply the function along - """ dim_indx: int = 0 def operate(self, data): - """ - Applies the operation for which the class is named for. + """Applies the operation for which the class is named for. Args: data: CSDM object - """ return data.fft(axis=self.dim_indx) @@ -118,5 +135,5 @@ class FFT(IFFT): pass -class complex_conjugate(abstractOperation): +class complex_conjugate(AbstractOperation): pass diff --git a/src/mrsimulator/signal_processing/_base.py b/src/mrsimulator/signal_processing/_base.py index 375d46ceb..d3dd21287 100644 --- a/src/mrsimulator/signal_processing/_base.py +++ b/src/mrsimulator/signal_processing/_base.py @@ -1,12 +1,14 @@ # -*- coding: utf-8 -*- -"""The abstractOperation class.""" +"""The AbstractOperation class.""" from mrsimulator.utils.parseable import Parseable -class abstractOperation(Parseable): - """ - A base class for signal processing operations - """ +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" + + +class AbstractOperation(Parseable): + """A base class for signal processing operations.""" @property def function(self): @@ -30,37 +32,3 @@ def parse_dict_with_units(cls, py_dict): if "function" in my_dict_copy.keys(): my_dict_copy.pop("function") return super().parse_dict_with_units(my_dict_copy) - - -# class FFT(abstractOperation): -# """ -# Class for applying a -# Fourier transform to a dependent variable of simulation data - -# dimension: int. Data dimension to apply the function along - -# """ - -# dimension: int = 0 - -# 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 data.fft(axis=self.dimension) - - -# class OperationList(BaseModel): -# """ -# Container class to hold a list of operations to be applied to a given -# dependent variable in the data - -# dependent_variable: The dependent variable of the data to apply operations to -# operations: List of abstractOperation based operations to sequentially apply to data. -# """ - -# dependent_variable: int = 0 -# operations: List[abstractOperation] diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index 2e9ec066e..d0a8023e5 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -1,6 +1,5 @@ # -*- coding: utf-8 -*- """The Event class.""" -from copy import deepcopy from sys import modules from typing import ClassVar from typing import Dict @@ -8,19 +7,19 @@ import numpy as np -from ._base import abstractOperation +from ._base import AbstractOperation +__author__ = "Maxwell C. Venetos" +__email__ = "maxvenetos@gmail.com" -class AbstractApodization(abstractOperation): + +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): - copy_dict = deepcopy(py_dict) - copy_dict.pop("type") - - obj = super().parse_dict_with_units(copy_dict) + obj = super().parse_dict_with_units(py_dict) return getattr(modules[__name__], py_dict["type"])(**obj.dict()) @property @@ -29,6 +28,7 @@ def function(self): @property def type(self): + """The type apodization function.""" return self.__class__.__name__ def set_property_units(self, unit, prop): @@ -80,9 +80,8 @@ def operate(self, data, fn, prop_name, prop_value): class Gaussian(AbstractApodization): - r""" - Class for applying a Gaussian function - to a dependent variable of simulation data + r"""Class for applying a Gaussian function to a dependent variable of simulation + data. .. math:: f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2} @@ -92,7 +91,6 @@ class Gaussian(AbstractApodization): sigma: float. Standard deviation of Gaussian function dep_var_indx: int. Data dependent variable index to apply the function to. If type None, will be applied to every dependent variable. - """ sigma: float = 0 @@ -119,9 +117,8 @@ def operate(self, data): class Exponential(AbstractApodization): - r""" - Class for applying an exponential - Lorentzian function to a dependent variable of simulation data + r"""Class for applying an exponential Lorentzian function to a dependent variable + of simulation data. .. math:: f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)} @@ -131,7 +128,6 @@ class Exponential(AbstractApodization): Lambda: float. Width parameter dep_var_indx: int. Data dependent variable index to apply the function to. If type None, will be applied to every dependent variable. - """ Lambda: float = 0 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..0966fb816 --- /dev/null +++ b/src/mrsimulator/signal_processing/tests/test_signal_processing.py @@ -0,0 +1,35 @@ +# -*- coding: utf-8 -*- +import mrsimulator.signal_processing as sp +import mrsimulator.signal_processing.apodization as apo + + +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 + + # 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 == post_sim_2 diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index d71d73dfd..4a9d8b028 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -55,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). diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 13ff24790..ba3b5f219 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -22,7 +22,7 @@ ["].shielding_symmetric.", "_shielding_symmetric_"], ["].quadrupolar.", "_quadrupolar_"], ["].abundance", "_abundance"], - ["methods[", "METHODS_"], # why does methods needs to be parameterized? + ["methods[", "METHODS_"], # why does methods need to be parameterized? ["].post_simulation", "_POST_SIM_"], [".scale", "scale"], [".apodization[", "APODIZATION_"], 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) """ From 0e7d4f0d2ab8a446cef992249f159ec231fae42c Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 6 Jul 2020 21:12:56 -0400 Subject: [PATCH 15/26] add and update docs for operations. --- conftest.py | 4 ++ .../operations/operation_documentation.rst | 15 ++++++ docs/api_py/operations/operations_summary.rst | 49 +++++++++++++++++++ docs/api_py/py_api.rst | 1 + ...y => plot_5_post_sim_signal_processing.py} | 5 +- src/mrsimulator/signal_processing/__init__.py | 40 ++++++++++++--- .../signal_processing/apodization.py | 49 +++++++++++++------ src/mrsimulator/spectral_fitting.py | 26 ++++------ 8 files changed, 149 insertions(+), 40 deletions(-) create mode 100644 docs/api_py/operations/operation_documentation.rst create mode 100644 docs/api_py/operations/operations_summary.rst rename examples/1D_simulation/{plot_zpost_sim_signal_processing.py => plot_5_post_sim_signal_processing.py} (98%) 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/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 8adf26a3e..cd76b8d43 100644 --- a/docs/api_py/py_api.rst +++ b/docs/api_py/py_api.rst @@ -16,6 +16,7 @@ Python-API References methods fitting signal_processing + operations/operations_summary .. parameterized_tensor .. interaction diff --git a/examples/1D_simulation/plot_zpost_sim_signal_processing.py b/examples/1D_simulation/plot_5_post_sim_signal_processing.py similarity index 98% rename from examples/1D_simulation/plot_zpost_sim_signal_processing.py rename to examples/1D_simulation/plot_5_post_sim_signal_processing.py index 04fd5b569..51f5435da 100644 --- a/examples/1D_simulation/plot_zpost_sim_signal_processing.py +++ b/examples/1D_simulation/plot_5_post_sim_signal_processing.py @@ -5,7 +5,6 @@ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. sectionauthor:: Maxwell C. Venetos """ -# sphinx_gallery_thumbnail_number = 5 # %% # After running a simulation, we often want to process the resulting spectrum. # For example, we may want to scale the intensities to match the experiment or @@ -13,6 +12,7 @@ # demonstrate the use of the `SignalProcessor` class to apply various operations # to simulation data. # global plot configuration +# sphinx_gallery_thumbnail_number = 5 import matplotlib as mpl import matplotlib.pyplot as plt @@ -23,8 +23,7 @@ # %% # We will create a hypothetical two-site Si simulation to illustrate post-simulation # signal processing. We will begin by processing the entire spectrum and follow up by -# decomposing the spectrum and processing each signal indepenently. - +# decomposing the spectrum and processing each signal independently. from mrsimulator import SpinSystem from mrsimulator import Simulator diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py index 1bebf0887..cf719f6fb 100644 --- a/src/mrsimulator/signal_processing/__init__.py +++ b/src/mrsimulator/signal_processing/__init__.py @@ -21,10 +21,10 @@ class SignalProcessor(BaseModel): ---------- data: CSDM object. - From simulation + The data from the simulation. operations: List - List of operation lists + A list of operations. Examples -------- @@ -93,7 +93,17 @@ def apply_operations(self, **kwargs): class Scale(AbstractOperation): """ - Class for applying a scaling factor to a dependent variable of simulation data. + 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 @@ -113,11 +123,16 @@ def operate(self, data): class IFFT(AbstractOperation): """ - Class for applying an inverse Fourier transform to a dependent variable of - simulation data. + Apply an inverse Fourier transform on the dependent variables from the simulation + data object. Args: - dim_indx: int. Data dimension to apply the function along + int dim_indx: Data dimension index along which the function is applied. + + Example + ------- + + >>> operation2 = sp.IFFT(dim_indx=0) """ dim_indx: int = 0 @@ -132,6 +147,19 @@ def operate(self, data): 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 diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index d0a8023e5..b9c41b5fc 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -80,17 +80,28 @@ def operate(self, data, fn, prop_name, prop_value): class Gaussian(AbstractApodization): - r"""Class for applying a Gaussian function to a dependent variable of simulation - data. + r"""Apodize a dependent variable of the simulation data object by a Gaussian + function. The function follows .. math:: - f(\vec{x}) = \vec{x}*e^{-2*(\vec{x} * \sigma * \pi)^2} + 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: - dim_indx: int. Data dimension to apply the function along. - sigma: float. Standard deviation of Gaussian function - dep_var_indx: int. Data dependent variable index to apply the function to. - If type None, will be applied to every dependent variable. + 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 @@ -101,7 +112,7 @@ class Gaussian(AbstractApodization): @staticmethod def fn(x, arg): - return np.exp(-((x * arg * np.pi) ** 2) * 2) + return np.exp(-2 * ((x * arg * np.pi) ** 2)) def operate(self, data): """ @@ -117,17 +128,25 @@ def operate(self, data): class Exponential(AbstractApodization): - r"""Class for applying an exponential Lorentzian function to a dependent variable - of simulation data. + r"""Apodize a dependent variable of the simulation data object by an exponential + function. The function follows .. math:: - f(\vec{x}) = \vec{x}*e^{-\Lambda * \abs{\vec{x}} * \pi)} + f(x) = e^{-\Lambda |{x}| \pi}, + + where :math:`x` are the coordinates of the data dimension and :math:`\Lambda` is + the width parameter. Args: - dim_indx: int. Data dimension to apply the function along. - Lambda: float. Width parameter - dep_var_indx: int. Data dependent variable index to apply the function to. - If type None, will be applied to every dependent variable. + int dim_indx: Data dimension index to apply the function along. + float Lambda: The width parameter, :math:`\Lambda`. + 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) """ Lambda: float = 0 diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index ba3b5f219..9e3550124 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -50,7 +50,6 @@ def _str_to_html(my_string): Returns: String object. - """ for item in ENCRYPTION_PAIRS: my_string = my_string.replace(*item) @@ -67,7 +66,6 @@ def _html_to_string(my_string): Returns: String Object. - """ for item in ENCRYPTION_PAIRS: my_string = my_string.replace(*item[::-1]) @@ -84,7 +82,6 @@ def _list_of_dictionaries(my_list): Returns: List Object. - """ return [item.dict() for item in my_list] @@ -96,13 +93,12 @@ 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): @@ -126,9 +122,11 @@ def _post_sim_LMFIT_params(post_sim): Creates an LMFIT Parameters object for SignalProcessor operations involved in spectrum fitting - post_sim: SignalProcessor object + Args: + post_sim: SignalProcessor object - returns: Parameters object + Returns: + Parameters object """ temp_dict = {} # for item in post_sim.operations: @@ -155,13 +153,11 @@ def _post_sim_LMFIT_params(post_sim): def _update_post_sim_from_LMFIT_params(params, post_sim): - """ - Updates SignalProcessor operation arguments from an - LMFIT Parameters object - - params: LMFIT Parameters object - post_sim: SignalProcessor object + """Updates SignalProcessor operation arguments from an LMFIT Parameters object + Args: + params: LMFIT Parameters object + post_sim: SignalProcessor object """ temp_dict = {} arg_dict = {"Gaussian": "sigma", "Exponential": "Lambda", "Scale": "factor"} @@ -204,7 +200,6 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): Returns: LMFIT Parameters object. - """ if not FOUND_LMFIT: error = ( @@ -281,7 +276,6 @@ def LMFIT_min_function(params, sim, post_sim=None): Returns: Array of the differences between the simulation and the experimental data. - """ if not isinstance(params, Parameters): raise ValueError( From dbcf568c076c066b331c68de851b4baa9c36c444 Mon Sep 17 00:00:00 2001 From: mVenetos97 <53448901+mVenetos97@users.noreply.github.com> Date: Tue, 7 Jul 2020 12:19:36 -0400 Subject: [PATCH 16/26] spectral fitting tests --- .../tests/test_spectral_fitting.py | 104 ++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 src/mrsimulator/signal_processing/tests/test_spectral_fitting.py 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..3717b1ad4 --- /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 = "ISO_0_SITES_0_isotropic_chemical_shift" + + str2 = "spin_systems[0].sites[0].quadrupolar.Cq" + str2_encoded = "ISO_0_SITES_0_quadrupolar_Cq" + + assert sf._str_to_html(str1) == str1_encoded + assert sf._str_to_html(str2) == str2_encoded + + +def test_02(): + # html_to_str + str1 = "spin_systems[0].sites[0].isotropic_chemical_shift" + str1_encoded = "ISO_0_SITES_0_isotropic_chemical_shift" + + str2 = "spin_systems[0].sites[0].quadrupolar.Cq" + str2_encoded = "ISO_0_SITES_0_quadrupolar_Cq" + + assert sf._html_to_string(str1_encoded) == str1 + assert sf._html_to_string(str2_encoded) == str2 + + +def test_03(): + # post_sim_LMFIT_params + op_list = [ + sp.IFFT(dim_indx=0), + apo.Exponential(Lambda=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(Lambda=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].Lambda == 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(Lambda=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 = { + "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, + "operation_1_Exponential": 100, + "operation_3_Scale": 10, + } + + assert params.valuesdict() == valuesdict, "Parameter creation failed" From 2676f740761a5b2983cb2e4ad648398d8ed376c3 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Tue, 7 Jul 2020 19:52:14 -0400 Subject: [PATCH 17/26] fixed lgtm issue. --- src/mrsimulator/signal_processing/apodization.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index b9c41b5fc..64e5c1c92 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -55,7 +55,7 @@ def _get_dv_indexes(indexes, n): if isinstance(indexes, (list, tuple)): return np.asarray(indexes) - def operate(self, data, fn, prop_name, prop_value): + def _operate(self, data, fn, prop_name, prop_value): """A generic operation function. Args: @@ -122,9 +122,7 @@ def operate(self, data): dep_var: int. The index of the dependent variable to apply operation to """ - return super().operate( - data, fn=self.fn, prop_name="sigma", prop_value=self.sigma - ) + return self._operate(data, fn=self.fn, prop_name="sigma", prop_value=self.sigma) class Exponential(AbstractApodization): @@ -167,6 +165,6 @@ def operate(self, data): dep_var: int. The index of the dependent variable to apply operation to """ - return super().operate( + return self._operate( data, fn=self.fn, prop_name="Lambda", prop_value=self.Lambda ) From fd27abbbc810ac197bc86b3f47998ac5455ba4c5 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Wed, 8 Jul 2020 12:48:02 -0400 Subject: [PATCH 18/26] enforce the max-line limit to 88. fixed HTML render for fitting examples. fixed thumbnail image for line-broadening example --- .flake8 | 4 +- .../plot_5_post_sim_signal_processing.py | 4 +- .../plot_1_mrsimFitExample_cuspidine.py | 81 ++++++++----------- .../Fitting/plot_2_mrsimFitExample_O17.py | 26 +++--- src/mrsimulator/signal_processing/__init__.py | 7 +- src/mrsimulator/spectral_fitting.py | 9 ++- 6 files changed, 64 insertions(+), 67 deletions(-) diff --git a/.flake8 b/.flake8 index d7840d4ad..696d05318 100644 --- a/.flake8 +++ b/.flake8 @@ -1,7 +1,9 @@ [flake8] -ignore = E402 #C901 E265 E402 E501 +ignore = E402 exclude = .eggs, *.egg,build, src/mrsimulator/__init__.py filename = *.pyx, *py max-line-length = 88 max-complexity = 12 select = C,E,F,W,N8 +count = True +statistics = True diff --git a/examples/1D_simulation/plot_5_post_sim_signal_processing.py b/examples/1D_simulation/plot_5_post_sim_signal_processing.py index 51f5435da..86f6ebf6b 100644 --- a/examples/1D_simulation/plot_5_post_sim_signal_processing.py +++ b/examples/1D_simulation/plot_5_post_sim_signal_processing.py @@ -11,14 +11,14 @@ # apodize the signal to simulate line broadening. The following example will # demonstrate the use of the `SignalProcessor` class to apply various operations # to simulation data. -# global plot configuration -# sphinx_gallery_thumbnail_number = 5 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.25, 3.0] +# sphinx_gallery_thumbnail_number = 4 # %% # We will create a hypothetical two-site Si simulation to illustrate post-simulation diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index d6608daac..988a6c071 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -5,19 +5,21 @@ ^^^^^^^^^^^^^^^^^^^ .. 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]_ +# 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]_ # # We will begin by importing *matplotlib* and establishing figure size. import matplotlib as mpl @@ -26,9 +28,11 @@ font = {"weight": "light", "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. +# Next we will import `csdmpy `_ +# and loading the data file. import csdmpy as cp filename = "https://osu.box.com/shared/static/a45xj96iekdjrs2beri0nkrow4vjewdh.csdf" @@ -42,15 +46,13 @@ 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: - +# 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 SpinSystem from mrsimulator import Simulator # %% # **Step 1** Create the site. - site = { "isotope": "29Si", "isotropic_chemical_shift": "-80.0 ppm", @@ -59,7 +61,6 @@ # %% # **Step 2** Create the spin system for the site. - spin_system = { "name": "Si Site", "description": "A 29Si site in cuspidine", @@ -70,7 +71,6 @@ # %% # **Step 3** Create the Bloch Decay method. - from mrsimulator.methods import BlochDecaySpectrum method = BlochDecaySpectrum( @@ -94,7 +94,6 @@ # %% # **Step 4** Create the Simulator object and add the method and spin-system objects. - sim = Simulator() sim.spin_systems += [system_object] sim.methods += [method] @@ -103,12 +102,10 @@ # %% # **Step 5** simulate the spectrum. - sim.run() # %% # **Step 6** Create a SignalProcessor - import mrsimulator.signal_processing as sp import mrsimulator.signal_processing.apodization as apo @@ -121,10 +118,8 @@ post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) - # %% # ** Step 7** Process and plot the spectrum. - processed_data = post_sim.apply_operations() ax = plt.subplot(projection="csdm") @@ -134,16 +129,13 @@ 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. - -# %% - - +# 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, fit_report params = Parameters() @@ -156,9 +148,10 @@ params # %% -# 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. +# 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. def test_function(params, sim, post_sim): @@ -184,21 +177,19 @@ def test_function(params, sim, post_sim): # %% -# 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 magic angle spinning fitting is we must use a discrete minimization method -# such as 'powell' as the chemical shift varies discretely +# 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 magic angle spinning fitting is we must use a discrete minimization +# method such as 'powell' as the chemical shift varies discretely # %% **Step 8** Perform minimization. - minner = Minimizer(test_function, params, fcn_args=(sim, post_sim)) result = minner.minimize(method="powell") print(fit_report(result)) # %% # **Step 9** Plot the fitted spectrum. - - plt.figsize = (4, 3) residual = synthetic_experiment.copy() residual[:] = result.residual @@ -215,12 +206,8 @@ def test_function(params, sim, post_sim): 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 489a1901a..4bb5f9f73 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -5,17 +5,18 @@ ^^^^^^^^^^^^^^^^^^^^^^^^ .. 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. +# 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]_ # # We will begin by importing *matplotlib* and establishing figure size. @@ -25,9 +26,11 @@ font = {"weight": "light", "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. +# Next we will import `csdmpy `_ +# and loading the data file. import csdmpy as cp @@ -142,7 +145,8 @@ # %% # 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*. +# fitting. We will be using the +# `Parameters `_ class from *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 diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py index cf719f6fb..520842180 100644 --- a/src/mrsimulator/signal_processing/__init__.py +++ b/src/mrsimulator/signal_processing/__init__.py @@ -29,7 +29,7 @@ class SignalProcessor(BaseModel): Examples -------- - >>> post_sim = SignalProcessor(data = csdm_object, operations = [op1, op2]) # doctest: +SKIP + >>> post_sim = SignalProcessor(data=csdm_data, operations=[o1, o2]) # doctest: +SKIP """ data: cp.CSDM = None @@ -62,8 +62,9 @@ def parse_dict_with_units(self, py_dict): 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. + 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. diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 9e3550124..8a3ea9382 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -162,7 +162,8 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): temp_dict = {} arg_dict = {"Gaussian": "sigma", "Exponential": "Lambda", "Scale": "factor"} for param in params: - # iterating through the parameter list looking for only DEP_VAR (ie post_sim 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) @@ -172,7 +173,8 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): # 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 + # 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 @@ -181,7 +183,8 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): # 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 + # 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]]}" From df8e1a0a3715258c49531d0d102a2d9c4ed04823 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Thu, 9 Jul 2020 22:45:46 -0400 Subject: [PATCH 19/26] Fixed bug in spectral_fitting. Expanded on the description of the fitting example. Expanded on the description of the signal processing example. Fixed the title Sodium silicate -> Crystalline Sodium Metasilicate --- .../plot_5_post_sim_signal_processing.py | 47 ++--- .../Fitting/plot_2_mrsimFitExample_O17.py | 165 ++++++++++-------- src/mrsimulator/spectral_fitting.py | 4 +- 3 files changed, 109 insertions(+), 107 deletions(-) diff --git a/examples/1D_simulation/plot_5_post_sim_signal_processing.py b/examples/1D_simulation/plot_5_post_sim_signal_processing.py index 86f6ebf6b..566db1d71 100644 --- a/examples/1D_simulation/plot_5_post_sim_signal_processing.py +++ b/examples/1D_simulation/plot_5_post_sim_signal_processing.py @@ -13,6 +13,11 @@ # to simulation data. 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 SpinSystem +from mrsimulator.methods import BlochDecaySpectrum # global plot configuration font = {"weight": "light", "size": 9} @@ -24,10 +29,7 @@ # We will create a hypothetical two-site Si simulation to illustrate post-simulation # signal processing. We will begin by processing the entire spectrum and follow up by # decomposing the spectrum and processing each signal independently. -from mrsimulator import SpinSystem -from mrsimulator import Simulator - -# %% +# # **Step 1** Create the sites and spin system site1 = { @@ -49,7 +51,6 @@ "abundance": "100%", } - spin_system_2 = { "name": "site B", "description": "A test 29Si site", @@ -62,15 +63,11 @@ # %% # **Step 2** Create the simulation and add the spin system objects. - sim = Simulator() sim.spin_systems += [system_object_1, system_object_2] # %% # **Step 3** Create a Bloch decay spectrum method. - -from mrsimulator.methods import BlochDecaySpectrum - method_dict = { "channels": ["29Si"], "magnetic_flux_density": "7.1 T", @@ -87,15 +84,12 @@ # The above method is set up to record the :math:`^{29}\text{Si}` resonances at the # magic angle, spinning at 780 Hz and 7.1 T external magnetic flux density. The # resonances are recorded over 25 kHz spectral width using 2048 points. - -# %% +# # **Step 4** Simulate the spectra. - sim.run() # %% # **Step 5** Plot the spectrum. - ax = plt.subplot(projection="csdm") ax.plot(sim.methods[0].simulation, color="black", linewidth=1) ax.invert_xaxis() @@ -103,18 +97,19 @@ plt.show() # %% -# **Step 6** Create Post Simulation. +# **Step 6** Create the post simulation operation by first creating a list of +# processing operations -import mrsimulator.signal_processing as sp -import mrsimulator.signal_processing.apodization as apo - -# create list of processing operations +# list of processing operations op_list1 = [ sp.IFFT(dim_indx=0), apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), sp.FFT(dim_indx=0), ] +# %% +# and then create the :class:`~mrsimulator.signal_processing.SignalProcessor` class +# object as follows post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list1) # %% @@ -125,15 +120,13 @@ # Next, a Gaussian exponential with a broadening factor of 100 seconds is # multipled to the data and then finally another Fourier transform is applied # to the data to convert the data back to frequency space. - -# %% +# +# # **Step 7** Apply the signal processing. - processed_data = post_sim.apply_operations() # %% # **Step 8** Plot the processed spectrum. - ax = plt.subplot(projection="csdm") ax.plot(processed_data, color="black", linewidth=1) ax.invert_xaxis() @@ -146,16 +139,13 @@ # seperated. 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() # %% # **Step 5** Plot spectrum. - x, y0, y1 = sim.methods[0].simulation.to_list() plt.plot(x, y0, x, y1) @@ -167,7 +157,6 @@ # %% # **Step 6** Create Post Simulation. - op_list2 = [ sp.IFFT(dim_indx=0), apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), @@ -185,12 +174,10 @@ # %% # **Step 7** Apply the singla processing. - processed_data = post_sim.apply_operations() # %% # **Step 8** Plot the processed spectrum - x, y0, y1 = processed_data.to_list() plt.plot(x, y0, x, y1) @@ -199,5 +186,3 @@ plt.grid(color="gray", linestyle="--", linewidth=0.5, alpha=0.5) plt.tight_layout() plt.show() - -# %% diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 4bb5f9f73..be8c46843 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -1,8 +1,8 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Fitting Sodium Silicate. -^^^^^^^^^^^^^^^^^^^^^^^^ +Fitting Crystalline Sodium Metasilicate +======================================= .. sectionauthor:: Maxwell C. Venetos """ # %% @@ -19,68 +19,87 @@ # obtained :math:`^{17}\text{O}` spectrum. # We use the :math:`^{17}\text{O}` tensor information from Grandinetti `et. al.` [#f5]_ # -# We will begin by importing *matplotlib* and establishing figure size. +# 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 -font = {"weight": "light", "size": 9} +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 +# ------------------ +# +# **Step 0** +# Import the experimental data. In this example, we will import the data file serialized +# with the CSDM file-format. We use the +# `csdmpy `_ module to load the +# dataset. +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") +# plot of the dataset. ax = plt.subplot(projection="csdm") ax.plot(oxygen_experiment, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() - # %% +# Create Simulator object +# ----------------------- +# # 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 SpinSystem -from mrsimulator import Simulator - -# %% -# **Step 1** Create the site. - -O17_1 = { - "isotope": "17O", - "isotropic_chemical_shift": "60.0 ppm", - "quadrupolar": {"Cq": "4200000 Hz", "eta": 0.5}, -} +# spectrum. We will start by creating a guess ``SpinSystem`` objects. +# +# **Step 1** Create the guess sites. +O17_1 = Site( + isotope="17O", + isotropic_chemical_shift=60.0, # in ppm, + quadrupolar={"Cq": 4.2e6, "eta": 0.5}, # Cq in Hz +) -O17_2 = { - "isotope": "17O", - "isotropic_chemical_shift": "40.0 ppm", - "quadrupolar": {"Cq": "2400000 Hz", "eta": 0.5}, -} +O17_2 = Site( + isotope="17O", + isotropic_chemical_shift=40.0, # in ppm, + quadrupolar={"Cq": 2.4e6, "eta": 0.5}, # Cq in Hz +) # %% -# **Step 2** Create the spin system for the site. - -spin_system = {"sites": [O17_1, O17_2], "abundance": "100%"} # from the above code - -system_object = SpinSystem.parse_dict_with_units(spin_system) +# **Step 2** Create the spin systems for the guessed sites. +system_object = [SpinSystem(sites=[s]) for s in [O17_1, O17_2]] # from the above code # %% -# **Step 3** Create the Bloch Decay method. - -from mrsimulator.methods import BlochDecayCentralTransitionSpectrum - +# **Step 3** Create the method object. Note, when fitting, 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 experiment dimension metadata. 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 @@ -99,81 +118,79 @@ ) # %% -# **Step 4** Create the Simulator object and add the method and spin-system objects. +# Now we add the experimental data to the ``experiment`` attribute of the above method. +method.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] - -sim.methods[0].experiment = oxygen_experiment +sim.spin_systems = system_object +sim.methods = [method] # %% -# **Step 5** simulate the spectrum. - +# **Step 5** Simulate the spectrum. for iso in sim.spin_systems: # To avoid querying at every iteration we will save the relevant transition pathways iso.transition_pathways = method.get_transition_pathways(iso).tolist() sim.run() # %% -# **Step 6** Create a SignalProcessor - -import mrsimulator.signal_processing as sp -import mrsimulator.signal_processing.apodization as apo - -factor = oxygen_experiment.dependent_variables[0].components[0].max().real / 4 - +# **Step 6** Create the SignalProcessor class object and add to it the list of +# operations that will be applied to the simulated spectrum. The generic list of +# operations in NMR includes the line-broadening function. In the following, we add a +# scaling and a Lorentzian line-broadening function. Here, the Lorentzian +# line-broadening is defined as an exponential apodization operation sandwiched +# between two Fourier transformations. +factor = oxygen_experiment.max().real op_list = [ - sp.IFFT(dim_indx=0), - apo.Exponential(Lambda=100, dim_indx=0, dep_var_indx=0), - sp.FFT(dim_indx=0), + sp.IFFT(), + apo.Exponential(Lambda=100), + sp.FFT(), sp.Scale(factor=factor), ] - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) # %% -# ** Step 7** Process and plot the spectrum. - +# **Step 7** Process and plot the spectrum. processed_data = post_sim.apply_operations() ax = plt.subplot(projection="csdm") -ax.plot(processed_data, color="black", linewidth=1) +ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() plt.tight_layout() plt.show() # %% +# Least-squares minimization with LMFIT +# ------------------------------------- +# # 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*. # # 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` - -# %% +# :func:`~mrsimulator.spectral_fitting.make_LMFIT_parameters` +# # **Step 8** Create a list of parameters to vary during fitting. - - -from mrsimulator.spectral_fitting import make_LMFIT_parameters - params = make_LMFIT_parameters(sim, post_sim) -params.pretty_print() -# %% -# **Step 9** Perform minimization. -from mrsimulator.spectral_fitting import LMFIT_min_function -from lmfit import Minimizer, report_fit +# %% +# **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["ISO_0_abundance"].vary = False # fix the abundance +params.pretty_print() +# %% +# **Step 9** Perform the minimization. minner = Minimizer(LMFIT_min_function, params, fcn_args=(sim, post_sim)) result = minner.minimize() report_fit(result) # %% # **Step 10** Plot the fitted spectrum. - - plt.figsize = (4, 3) residual = oxygen_experiment.copy() residual[:] = result.residual diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 8a3ea9382..9f454dd4d 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -223,8 +223,8 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): # expression for the last abundance. last_abund = f"{length - 1}_abundance" - expression = "100" + "-".join([f"{START}{i}_abundance" for i in range(length - 1)]) - + 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: params.add( From 106f37a2c435277917d25163e92e7ce61744d487 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Tue, 21 Jul 2020 18:44:26 -0400 Subject: [PATCH 20/26] removed signal_processing py file added signal_processing rst docs. Update example files. --- docs/index.rst | 1 + docs/signal_processing.rst | 281 ++++++++++++++++++ .../plot_5_post_sim_signal_processing.py | 188 ------------ .../plot_1_mrsimFitExample_cuspidine.py | 82 +++-- .../Fitting/plot_2_mrsimFitExample_O17.py | 2 +- 5 files changed, 323 insertions(+), 231 deletions(-) create mode 100644 docs/signal_processing.rst delete mode 100644 examples/1D_simulation/plot_5_post_sim_signal_processing.py diff --git a/docs/index.rst b/docs/index.rst index 2d8428f83..8fa01e849 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/signal_processing.rst b/docs/signal_processing.rst new file mode 100644 index 000000000..cb1bd33cc --- /dev/null +++ b/docs/signal_processing.rst @@ -0,0 +1,281 @@ + +Post Simulation Signal Processing +================================= +.. sectionauthor:: Maxwell C. Venetos + + +After running a simulation, we often want to process the resulting spectrum. +For example, we may want to scale the intensities to match the experiment or +apodize the signal to simulate line-broadening. The following section will +demonstrate the use of the `SignalProcessor` class in applying various operations +to the simulation data. + +Import the relevent modules. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> import matplotlib as mpl + >>> import matplotlib.pyplot as plt + >>> import mrsimulator.signal_processing as sp + >>> import mrsimulator.signal_processing.apodization as apo + ... + >>> # 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:`previous section ` 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 post-simulation signal processing steps. + +.. 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": -5000, # 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 +-------------------------- + +Processing single spectrum +'''''''''''''''''''''''''' + +.. We will begin by illustrating the signal processing of a entire spectrum, that is, +.. the combined spectrum arsing from the above two spin-systems. + +All signal processing operations are located in the `signal_processing` module of the +``mrsimulator`` library. Within the module is the `apodization` sub-module, which is +used in apodizing the signal. An apodization is a point-wise multiplication operation of +the input signal with the apodization method. Please read our :ref:`operations_api` +documentation for a complete list of operations. + +Let's see how we can use this module and its 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 + +The signal processing operations are given as an ordered list of operations. In this +workflow, the result from the previous operation becomes the input for the next operation. +The following example shows an application of a single operation---a Gaussian +line-broadening. Note, for almost all NMR spectrum simulation, the +post-simulation processing is a convolution, including the line-broadening. You may +achieve this convolution by sandwiching the corresponding apodization operation between +the two Fourier transformations. In this case, the Gaussian convolution is defined as +a set of three operations, as follows, + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> # list of processing operations + >>> op_list1 = [sp.IFFT(), apo.Gaussian(sigma=100), sp.FFT()] + +The above signal processing procedure is set to sequentially apply the list of +operations to simulate a Gaussian broadening of the spectrum. +The operation procedure will first perform an inverse Fourier Transform to convert the +frequency domain data 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 resulting time-domain signal back to the frequency domain. + +To apply the above list of operations to the simulation data, use the +:class:`~mrsimulator.signal_processing.SignalProcessor` class instance as follows + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list1) + >>> processed_data = post_sim.apply_operations() + +Here, the first line creates an instance of the `SignalProcessor` class whose attributes, +`data` and `operations`, hold the simulation data and the list of operations, respectively. +The second line applies the operations and returns a 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 + + +Processing individual 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 are 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 + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP + +Because the simulation is stored as a CSDM object, each sub-spectrum is a +dependent-variable of the simulation dataset, each 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: + + >>> op_list2 = [ + ... sp.IFFT(), + ... apo.Gaussian(sigma=50, dep_var_indx=0), + ... apo.Exponential(Lambda=200, dep_var_indx=1), + ... sp.FFT(), + ... ] + +The above signal processing procedure will first apply 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, create the ``SignalProcessor`` class, add the data and the list of operations, +and apply the operations. + +.. plot:: + :format: doctest + :context: close-figs + :include-source: + + >>> post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list2) + >>> processed_data = post_sim.apply_operations() + +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 + >>> plt.tight_layout() # doctest: +SKIP + >>> plt.show() # doctest: +SKIP diff --git a/examples/1D_simulation/plot_5_post_sim_signal_processing.py b/examples/1D_simulation/plot_5_post_sim_signal_processing.py deleted file mode 100644 index 566db1d71..000000000 --- a/examples/1D_simulation/plot_5_post_sim_signal_processing.py +++ /dev/null @@ -1,188 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -""" -Post Simulation Signal Processing. -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -.. sectionauthor:: Maxwell C. Venetos -""" -# %% -# After running a simulation, we often want to process the resulting spectrum. -# For example, we may want to scale the intensities to match the experiment or -# apodize the signal to simulate line broadening. The following example will -# demonstrate the use of the `SignalProcessor` class to apply various operations -# to simulation data. -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 SpinSystem -from mrsimulator.methods import BlochDecaySpectrum - -# global plot configuration -font = {"weight": "light", "size": 9} -mpl.rc("font", **font) -mpl.rcParams["figure.figsize"] = [4.25, 3.0] -# sphinx_gallery_thumbnail_number = 4 - -# %% -# We will create a hypothetical two-site Si simulation to illustrate post-simulation -# signal processing. We will begin by processing the entire spectrum and follow up by -# decomposing the spectrum and processing each signal independently. -# -# **Step 1** Create the sites and spin system - -site1 = { - "isotope": "29Si", - "isotropic_chemical_shift": "-75.0 ppm", - "shielding_symmetric": {"zeta": "-60 ppm", "eta": 0.6}, -} - -site2 = { - "isotope": "29Si", - "isotropic_chemical_shift": "-80.0 ppm", - "shielding_symmetric": {"zeta": "-70 ppm", "eta": 0.5}, -} - -spin_system_1 = { - "name": "site A", - "description": "A test 29Si site", - "sites": [site1], # from the above code - "abundance": "100%", -} - -spin_system_2 = { - "name": "site B", - "description": "A test 29Si site", - "sites": [site2], # from the above code - "abundance": "100%", -} - -system_object_1 = SpinSystem.parse_dict_with_units(spin_system_1) -system_object_2 = SpinSystem.parse_dict_with_units(spin_system_2) - -# %% -# **Step 2** Create the simulation and add the spin system objects. -sim = Simulator() -sim.spin_systems += [system_object_1, system_object_2] - -# %% -# **Step 3** Create a Bloch decay spectrum method. -method_dict = { - "channels": ["29Si"], - "magnetic_flux_density": "7.1 T", - "rotor_angle": "54.735 deg", - "rotor_frequency": "780 Hz", - "spectral_dimensions": [ - {"count": 2048, "spectral_width": "25 kHz", "reference_offset": "-5 kHz"} - ], -} -method_object = BlochDecaySpectrum.parse_dict_with_units(method_dict) -sim.methods += [method_object] - -# %% -# The above method is set up to record the :math:`^{29}\text{Si}` resonances at the -# magic angle, spinning at 780 Hz and 7.1 T external magnetic flux density. The -# resonances are recorded over 25 kHz spectral width using 2048 points. -# -# **Step 4** Simulate the spectra. -sim.run() - -# %% -# **Step 5** Plot the spectrum. -ax = plt.subplot(projection="csdm") -ax.plot(sim.methods[0].simulation, color="black", linewidth=1) -ax.invert_xaxis() -plt.tight_layout() -plt.show() - -# %% -# **Step 6** Create the post simulation operation by first creating a list of -# processing operations - -# list of processing operations -op_list1 = [ - sp.IFFT(dim_indx=0), - apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), - sp.FFT(dim_indx=0), -] - -# %% -# and then create the :class:`~mrsimulator.signal_processing.SignalProcessor` class -# object as follows -post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list1) - -# %% -# The above signal processing procedure will apply the list of operations -# sequentially to the 0-index dependent variable to result in Gaussian -# broadening of the spectrum. The operation procedure will first perform -# a Fourier Transform to convert the data in frequency space to time space. -# Next, a Gaussian exponential with a broadening factor of 100 seconds is -# multipled to the data and then finally another Fourier transform is applied -# to the data to convert the data back to frequency space. -# -# -# **Step 7** Apply the signal processing. -processed_data = post_sim.apply_operations() - -# %% -# **Step 8** Plot the processed spectrum. -ax = plt.subplot(projection="csdm") -ax.plot(processed_data, color="black", linewidth=1) -ax.invert_xaxis() -plt.tight_layout() -plt.show() - -# %% -# The above code resulted in the same processing to be applied -# to both signals because in the simulation the signals were not -# seperated. 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() - -# %% -# **Step 5** Plot spectrum. -x, y0, y1 = sim.methods[0].simulation.to_list() - -plt.plot(x, y0, x, y1) -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) -plt.tight_layout() -plt.show() - -# %% -# **Step 6** Create Post Simulation. -op_list2 = [ - sp.IFFT(dim_indx=0), - apo.Gaussian(sigma=100, dim_indx=0, dep_var_indx=0), - apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=1), - sp.FFT(dim_indx=0), -] - -post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list2) - -# %% -# The above signal processing procedure will apply Gaussian -# broadening to dependent variable index 0 (-75 ppm Si site) -# and apply Lorentzian broadening to dependent variable -# index 1 (-80 ppm Si site). - -# %% -# **Step 7** Apply the singla processing. -processed_data = post_sim.apply_operations() - -# %% -# **Step 8** Plot the processed spectrum -x, y0, y1 = processed_data.to_list() - -plt.plot(x, y0, x, y1) -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) -plt.tight_layout() -plt.show() diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 988a6c071..90dde653e 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -1,8 +1,8 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Fitting Cusipidine. -^^^^^^^^^^^^^^^^^^^ +Fitting Cusipidine +================== .. sectionauthor:: Maxwell C. Venetos """ # %% @@ -25,12 +25,14 @@ import matplotlib as mpl import matplotlib.pyplot as plt -font = {"weight": "light", "size": 9} +font = {"size": 9} mpl.rc("font", **font) -mpl.rcParams["figure.figsize"] = [4.25, 3.0] +mpl.rcParams["figure.figsize"] = [4.5, 3.0] # sphinx_gallery_thumbnail_number = 3 # %% +# Import the dataset +# ------------------ # Next we will import `csdmpy `_ # and loading the data file. import csdmpy as cp @@ -46,33 +48,36 @@ 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 SpinSystem -from mrsimulator import Simulator +# Create a "fitting model" +# ------------------------ +# Before you can fit a simulation to the data, 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 # %% -# **Step 1** Create the site. -site = { - "isotope": "29Si", - "isotropic_chemical_shift": "-80.0 ppm", - "shielding_symmetric": {"zeta": "-60 ppm", "eta": 0.6}, -} +# **Step 1:** Create the guess sites and spin systems. Since it is most likely that you +# will be fitting for the spin-system parameters using some non-linear fitting +# algorithm, as a general recommendation, the guess spin-system(s) should be a good +# starting point. +site = dict( + isotope="29Si", + isotropic_chemical_shift=-80.0, # in ppm, + shielding_symmetric={"zeta": -60, "eta": 0.6}, # zeta in ppm +) -# %% -# **Step 2** Create the spin system for the site. -spin_system = { - "name": "Si Site", - "description": "A 29Si site in cuspidine", - "sites": [site], # from the above code - "abundance": "100%", -} -system_object = SpinSystem.parse_dict_with_units(spin_system) +system_object = SpinSystem( + name="Si Site", + description="A 29Si site in cuspidine", + sites=[site], # from the above code + abundance=100, +) # %% -# **Step 3** Create the Bloch Decay method. -from mrsimulator.methods import BlochDecaySpectrum - +# **Step 2:** Create the method. It is highly likely that the method used in the +# simulator and its parameters are well known. When creating the method object, set the +# value of the method parameters to the respective values used in the experiment. method = BlochDecaySpectrum( channels=["29Si"], magnetic_flux_density=7.1, # in T @@ -87,16 +92,10 @@ ) # %% -# The above method is set up to record the :math:`^{29}\text{Si}` resonances at the -# magic angle, spinning at 780 Hz and 7.1 T external magnetic flux density. -# The resonances are recorded over 25 kHz spectral width ofset by -5000 Hz and -# using 2046 points. - -# %% -# **Step 4** Create the Simulator object and add the method and spin-system objects. +# **Step 3:** Create the Simulator object and add the method and spin-system objects. sim = Simulator() -sim.spin_systems += [system_object] -sim.methods += [method] +sim.spin_systems = [system_object] +sim.methods = [method] sim.methods[0].experiment = synthetic_experiment @@ -105,23 +104,22 @@ sim.run() # %% -# **Step 6** Create a SignalProcessor +# **Step 6** Create a SignalProcessor and apply post simulation processing. import mrsimulator.signal_processing as sp import mrsimulator.signal_processing.apodization as apo op_list = [ - sp.IFFT(dim_indx=0), - apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=0), - sp.FFT(dim_indx=0), + sp.IFFT(), + apo.Exponential(Lambda=200), + sp.FFT(), sp.Scale(factor=1), ] post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) - -# %% -# ** Step 7** Process and plot the spectrum. processed_data = post_sim.apply_operations() +# %% +# **Step 7** The plot the spectrum. ax = plt.subplot(projection="csdm") ax.plot(processed_data, color="black", linewidth=1) ax.invert_xaxis() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index be8c46843..31f3210b3 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -141,7 +141,7 @@ # scaling and a Lorentzian line-broadening function. Here, the Lorentzian # line-broadening is defined as an exponential apodization operation sandwiched # between two Fourier transformations. -factor = oxygen_experiment.max().real +factor = oxygen_experiment.max() op_list = [ sp.IFFT(), apo.Exponential(Lambda=100), From bdf90c9f97833c3bb18d301deb4d57f71ac80688 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Thu, 23 Jul 2020 16:03:21 -0400 Subject: [PATCH 21/26] removed `data` attribute from SignalProcessor add `processed_data` attribute to SignalProcessor update docs relating to signal-processing update examples relating to signal processing fix SignalProcessor serialization bug showing after applying operation. --- docs/signal_processing.rst | 163 +++++++++++------ examples/Fitting/README.rst | 4 +- .../plot_1_mrsimFitExample_cuspidine.py | 166 +++++++++++------- .../Fitting/plot_2_mrsimFitExample_O17.py | 12 +- src/mrsimulator/signal_processing/__init__.py | 20 +-- .../signal_processing/apodization.py | 4 +- .../tests/test_apodization.py | 16 +- .../tests/test_signal_processing.py | 11 +- src/mrsimulator/spectral_fitting.py | 3 +- 9 files changed, 246 insertions(+), 153 deletions(-) diff --git a/docs/signal_processing.rst b/docs/signal_processing.rst index cb1bd33cc..5c2da3864 100644 --- a/docs/signal_processing.rst +++ b/docs/signal_processing.rst @@ -4,13 +4,21 @@ Post Simulation Signal Processing .. sectionauthor:: Maxwell C. Venetos -After running a simulation, we often want to process the resulting spectrum. -For example, we may want to scale the intensities to match the experiment or -apodize the signal to simulate line-broadening. The following section will -demonstrate the use of the `SignalProcessor` class in applying various operations -to the simulation data. - -Import the relevent modules. +After running a simulation, you may find yourself in need of some +post-simulation signal processing operations. +For example, you may want to scale the intensities to match the experiment or +add line-broadening. 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 @@ -32,7 +40,8 @@ Simulating spectrum ------------------- Please refer to the :ref:`previous section ` 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 post-simulation signal processing steps. +single-site spin-systems to illustrate the use of the post-simulation signal processing +module. .. plot:: :format: doctest @@ -98,17 +107,14 @@ The plot the spectrum is shown below. Post-simulating processing -------------------------- -Processing single spectrum -'''''''''''''''''''''''''' - -.. We will begin by illustrating the signal processing of a entire spectrum, that is, -.. the combined spectrum arsing from the above two spin-systems. +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, which is -used in apodizing the signal. An apodization is a point-wise multiplication operation of -the input signal with the apodization method. Please read our :ref:`operations_api` -documentation for a complete list of operations. +``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. Let's see how we can use this module and its operations. Import the module and sub-module as @@ -121,14 +127,11 @@ sub-module as >>> import mrsimulator.signal_processing as sp >>> import mrsimulator.signal_processing.apodization as apo -The signal processing operations are given as an ordered list of operations. In this -workflow, the result from the previous operation becomes the input for the next operation. -The following example shows an application of a single operation---a Gaussian -line-broadening. Note, for almost all NMR spectrum simulation, the -post-simulation processing is a convolution, including the line-broadening. You may -achieve this convolution by sandwiching the corresponding apodization operation between -the two Fourier transformations. In this case, the Gaussian convolution is defined as -a set of three operations, as follows, +The signal processing is a series of operations 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. In the following +example, we show the application of a single operation—-a Gaussian line-broadening--using +the :class:`~mrsimulator.signal_processing.SignalProcessor` class, as follows, .. plot:: :format: doctest @@ -136,30 +139,43 @@ a set of three operations, as follows, :include-source: >>> # list of processing operations - >>> op_list1 = [sp.IFFT(), apo.Gaussian(sigma=100), sp.FFT()] + >>> post_sim = sp.SignalProcessor( + ... operations=[ + ... sp.IFFT(), apo.Gaussian(sigma=100), sp.FFT() + ... ] + ... ) -The above signal processing procedure is set to sequentially apply the list of -operations to simulate a Gaussian broadening of the spectrum. -The operation procedure will first perform an inverse Fourier Transform to convert the -frequency domain data to the time domain. Next, the time domain signal is apodized by a +The required attribute of the ``SignalProcessor`` class, `operations`, holds the list of +operations that gets applied to the simulation dataset. In the above example, 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 resulting time-domain signal back to the frequency domain. +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 data, use the -:class:`~mrsimulator.signal_processing.SignalProcessor` class instance as follows +:meth:`~mrsimulator.signal_processing.SignalProcessor.apply_operations` method of the +instance as follows .. plot:: :format: doctest :context: close-figs :include-source: - >>> post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list1) - >>> processed_data = post_sim.apply_operations() + >>> processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) -Here, the first line creates an instance of the `SignalProcessor` class whose attributes, -`data` and `operations`, hold the simulation data and the list of operations, respectively. -The second line applies the operations and returns a processed data. The plot of the -processed signal is shown below. +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 of the +signal processing operations. The plot of the processed signal is shown below. .. plot:: :format: doctest @@ -173,8 +189,8 @@ processed signal is shown below. >>> plt.show() # doctest: +SKIP -Processing individual sub-spectra -''''''''''''''''''''''''''''''''' +Applying operation to the sub-spectra +''''''''''''''''''''''''''''''''''''' .. spectrum and follow up by decomposing the spectrum and processing each signal .. independently. @@ -230,8 +246,8 @@ The plot of the spectra are shown below >>> plt.tight_layout() # doctest: +SKIP >>> plt.show() # doctest: +SKIP -Because the simulation is stored as a CSDM object, each sub-spectrum is a -dependent-variable of the simulation dataset, each sharing the same frequency dimension. +Because the simulation is stored as a CSDM [#f1]_ object, each sub-spectrum is a +dependent-variable of the simulation dataset, 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 @@ -244,29 +260,30 @@ operations. :context: close-figs :include-source: - >>> op_list2 = [ - ... sp.IFFT(), - ... apo.Gaussian(sigma=50, dep_var_indx=0), - ... apo.Exponential(Lambda=200, dep_var_indx=1), - ... sp.FFT(), - ... ] + >>> post_sim = sp.SignalProcessor( + ... operations=[ + ... sp.IFFT(), # convert to time-domain + ... apo.Gaussian(sigma=50, dep_var_indx=0), + ... apo.Exponential(Lambda=200, dep_var_indx=1), + ... sp.FFT(), # convert to frequency-domain + ... ] + ... ) -The above signal processing procedure will first apply 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. +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, create the ``SignalProcessor`` class, add the data and the list of operations, -and apply the operations. +As before, apply the operations with the +:meth:`~mrsimulator.signal_processing.SignalProcessor.apply_operations` method. .. plot:: :format: doctest :context: close-figs :include-source: - >>> post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list2) - >>> processed_data = post_sim.apply_operations() + >>> processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) The plot of the processed spectrum is shown below. @@ -279,3 +296,35 @@ The plot of the processed spectrum is shown below. >>> ax.plot(processed_data, alpha=0.75) # doctest: +SKIP >>> plt.tight_layout() # doctest: +SKIP >>> plt.show() # doctest: +SKIP + + +Serializing the operations list +------------------------------- + +You may 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'}, + {'Lambda': '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/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 90dde653e..58574c3ee 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -6,20 +6,19 @@ .. sectionauthor:: Maxwell C. Venetos """ # %% -# 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 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 model parameters. +# +# 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 as mpl @@ -33,14 +32,17 @@ # %% # Import the dataset # ------------------ -# Next we will import `csdmpy `_ -# and loading the data file. +# 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") +# plot of the synthetic dataset. ax = plt.subplot(projection="csdm") ax.plot(synthetic_experiment, color="black", linewidth=1) ax.invert_xaxis() @@ -50,21 +52,26 @@ # %% # Create a "fitting model" # ------------------------ -# Before you can fit a simulation to the data, 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. +# +# 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 # %% -# **Step 1:** Create the guess sites and spin systems. Since it is most likely that you -# will be fitting for the spin-system parameters using some non-linear fitting -# algorithm, as a general recommendation, the guess spin-system(s) should be a good -# starting point. +# **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, # in ppm, - shielding_symmetric={"zeta": -60, "eta": 0.6}, # zeta in ppm + isotropic_chemical_shift=-82.0, # in ppm, + shielding_symmetric={"zeta": -63, "eta": 0.4}, # zeta in ppm ) system_object = SpinSystem( @@ -75,9 +82,10 @@ ) # %% -# **Step 2:** Create the method. It is highly likely that the method used in the -# simulator and its parameters are well known. When creating the method object, set the -# value of the method parameters to the respective values used in the experiment. +# **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, # in T @@ -104,7 +112,7 @@ sim.run() # %% -# **Step 6** Create a SignalProcessor and apply post simulation processing. +# **Step 6** Create a SignalProcessor class and apply post simulation processing. import mrsimulator.signal_processing as sp import mrsimulator.signal_processing.apodization as apo @@ -115,46 +123,63 @@ sp.Scale(factor=1), ] -post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) -processed_data = post_sim.apply_operations() +post_sim = sp.SignalProcessor(operations=op_list) +processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) # %% -# **Step 7** The plot the spectrum. +# **Step 7** The plot the spectrum. We also plot the synthetic dataset for comparison. ax = plt.subplot(projection="csdm") -ax.plot(processed_data, color="black", linewidth=1) +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.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. +# 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 LMFIT on its +# `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. , as follows 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="Lambda", value=200) -params.add(name="factor", 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="Lambda", value=post_sim.operations[1].Lambda) +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, sim, post_sim): +def minimization_function(params, sim, post_sim): values = params.valuesdict() + # 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 @@ -163,37 +188,46 @@ def test_function(params, sim, post_sim): site.shielding_symmetric.eta = values["eta"] site.shielding_symmetric.zeta = values["zeta"] - # here we run the simulation + # run the simulation sim.run() + # update the SignalProcessor parameter and apply line broadening. post_sim.operations[3].factor = values["factor"] post_sim.operations[1].Lambda = values["Lambda"] - # here we apodize the signal to simulate line broadening - processed_data = post_sim.apply_operations() + processed_data = post_sim.apply_operations(sim.methods[0].simulation) + # return the difference vector. return intensity - processed_data.dependent_variables[0].components[0].real # %% -# With the synthetic data, simulation, and the parameters we are ready to perform the +# .. note:: +# To automate the fitting process, we provide a function to automatically 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 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 magic angle spinning fitting is we must use a discrete minimization -# method such as 'powell' as the chemical shift varies discretely - -# %% **Step 8** Perform minimization. -minner = Minimizer(test_function, params, fcn_args=(sim, post_sim)) +# 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") print(fit_report(result)) # %% -# **Step 9** Plot the fitted spectrum. +# **The plot the fitted spectrum** 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.xlabel("Frequency / Hz") plt.gca().invert_xaxis() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 31f3210b3..95db62e22 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -65,8 +65,8 @@ plt.show() # %% -# Create Simulator object -# ----------------------- +# Create a "fitting model" +# ------------------------ # # Next, we will want to create a ``simulator`` object that we will use to fit to our # spectrum. We will start by creating a guess ``SpinSystem`` objects. @@ -85,7 +85,7 @@ ) # %% -# **Step 2** Create the spin systems for the guessed sites. +# **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 # %% @@ -96,7 +96,7 @@ # 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 experiment dimension metadata. The remaining attribute values are set to +# the dimension metadata of the dimension. The remaining attribute values are set to # the experimental conditions. # get the count, increment, and coordinates_offset info from the experiment dimension. @@ -148,11 +148,11 @@ sp.FFT(), sp.Scale(factor=factor), ] -post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=op_list) +post_sim = sp.SignalProcessor(operations=op_list) # %% # **Step 7** Process and plot the spectrum. -processed_data = post_sim.apply_operations() +processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) ax = plt.subplot(projection="csdm") ax.plot(processed_data.real, color="black", linewidth=1) diff --git a/src/mrsimulator/signal_processing/__init__.py b/src/mrsimulator/signal_processing/__init__.py index 520842180..df3250531 100644 --- a/src/mrsimulator/signal_processing/__init__.py +++ b/src/mrsimulator/signal_processing/__init__.py @@ -14,25 +14,22 @@ class SignalProcessor(BaseModel): """ - Signal processing class to apply lists of various operations to individual - dependent variables of the data. + Signal processing class to apply a series of operations to the dependent variables + of the simulation dataset. Attributes ---------- - data: CSDM object. - The data from the simulation. - operations: List A list of operations. Examples -------- - >>> post_sim = SignalProcessor(data=csdm_data, operations=[o1, o2]) # doctest: +SKIP + >>> post_sim = SignalProcessor(operations=[o1, o2]) # doctest: +SKIP """ - data: cp.CSDM = None + processed_data: cp.CSDM = None operations: List[AbstractOperation] = [] class Config: @@ -63,7 +60,7 @@ def parse_dict_with_units(self, py_dict): 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 + object, where physical quantities are represented as string with a value and a unit. Returns: @@ -77,7 +74,7 @@ def to_dict_with_units(self): op["operations"] = lst return op - def apply_operations(self, **kwargs): + 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. @@ -85,9 +82,12 @@ def apply_operations(self, **kwargs): Returns: CSDM object: A copy of the data member with the operations applied to it. """ - copy_data = self.data.copy() + 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 diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index 64e5c1c92..65ea98a4d 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -36,9 +36,9 @@ 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 = dict(prop="s") + self.property_units = {f"{prop}": "s"} return - self.property_units = dict(prop="Hz") + self.property_units = {f"{prop}": "Hz"} @staticmethod def _get_dv_indexes(indexes, n): diff --git a/src/mrsimulator/signal_processing/tests/test_apodization.py b/src/mrsimulator/signal_processing/tests/test_apodization.py index 78bf9ac99..39d1e79d2 100644 --- a/src/mrsimulator/signal_processing/tests/test_apodization.py +++ b/src/mrsimulator/signal_processing/tests/test_apodization.py @@ -52,8 +52,8 @@ def test_scale(): - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_0) - data = post_sim.apply_operations() + 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() @@ -63,8 +63,8 @@ def test_scale(): def test_Lorentzian(): - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_1) - data = post_sim.apply_operations() + 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 @@ -78,8 +78,8 @@ def test_Lorentzian(): def test_Gaussian(): - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_2) - data = post_sim.apply_operations() + post_sim = sp.SignalProcessor(operations=PS_2) + data = post_sim.apply_operations(data=sim.methods[0].simulation) _, y0, y1, _ = data.to_list() sigma = 20 @@ -91,8 +91,8 @@ def test_Gaussian(): ), "Gaussian apodization amplitude failed" # test None for dep_var_indx - post_sim = sp.SignalProcessor(data=sim.methods[0].simulation, operations=PS_3) - data = post_sim.apply_operations() + 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 diff --git a/src/mrsimulator/signal_processing/tests/test_signal_processing.py b/src/mrsimulator/signal_processing/tests/test_signal_processing.py index 0966fb816..9308fc31c 100644 --- a/src/mrsimulator/signal_processing/tests/test_signal_processing.py +++ b/src/mrsimulator/signal_processing/tests/test_signal_processing.py @@ -1,6 +1,9 @@ # -*- 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(): @@ -13,6 +16,12 @@ def test_01(): 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_ == { @@ -32,4 +41,4 @@ def test_01(): # parse dict with units post_sim_2 = sp.SignalProcessor.parse_dict_with_units(dict_) - assert post_sim == post_sim_2 + assert post_sim.operations == post_sim_2.operations diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index 9f454dd4d..a105ddbdd 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -302,8 +302,7 @@ def LMFIT_min_function(params, sim, post_sim=None): _update_post_sim_from_LMFIT_params(params, post_sim) sim.run() - post_sim.data = sim.methods[0].simulation - processed_data = post_sim.apply_operations() + processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) # residual = np.asarray([]) if sim.config.decompose_spectrum == "spin_system": From a96c3e15f7e7b5699a673a405d45343d884c8f86 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sat, 25 Jul 2020 17:48:29 -0400 Subject: [PATCH 22/26] update docs for signal processing update examples to include post-simulation processing. changed ENCRYPTION_PAIRS to ENCODING_PAIRS changed _str_to_html to _str_encode changed _html_to_string to _str_decode changed Lambda to FWHM for exp apodization --- docs/signal_processing.rst | 67 +++++++++---------- docs/using_mrsimulator_objects.rst | 3 - examples/1D_simulation/plot_0_Wollastonite.py | 13 +++- .../1D_simulation/plot_1_PotassiumSulfate.py | 13 +++- examples/1D_simulation/plot_2_Coesite.py | 13 +++- .../plot_1_mrsimFitExample_cuspidine.py | 36 +++++----- .../Fitting/plot_2_mrsimFitExample_O17.py | 2 +- .../signal_processing/apodization.py | 63 ++++++++++++++--- .../tests/test_apodization.py | 10 +-- .../tests/test_spectral_fitting.py | 16 ++--- src/mrsimulator/spectral_fitting.py | 30 ++++----- 11 files changed, 162 insertions(+), 104 deletions(-) diff --git a/docs/signal_processing.rst b/docs/signal_processing.rst index 5c2da3864..33dd04c76 100644 --- a/docs/signal_processing.rst +++ b/docs/signal_processing.rst @@ -4,15 +4,14 @@ Post Simulation Signal Processing .. sectionauthor:: Maxwell C. Venetos -After running a simulation, you may find yourself in need of some -post-simulation signal processing operations. -For example, you may want to scale the intensities to match the experiment or -add line-broadening. 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. +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 @@ -27,8 +26,6 @@ operations to the simulation data. >>> import matplotlib as mpl >>> import matplotlib.pyplot as plt - >>> import mrsimulator.signal_processing as sp - >>> import mrsimulator.signal_processing.apodization as apo ... >>> # global plot configuration >>> font = {"weight": "light", "size": 9} @@ -38,7 +35,7 @@ operations to the simulation data. Simulating spectrum ------------------- -Please refer to the :ref:`previous section ` for a detailed description +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. @@ -77,7 +74,7 @@ module. ... { ... "count": 2048, ... "spectral_width": 25000, # in Hz - ... "reference_offset": -5000, # in Hz + ... "reference_offset": -5070, # in Hz ... "label": "$^{29}$Si resonances", ... } ... ], @@ -107,6 +104,10 @@ The plot the spectrum is shown below. 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 '''''''''''''''''''''''''''' @@ -116,8 +117,7 @@ apodization is a point-wise multiplication operation of the input signal with th apodizing vector. Please read our :ref:`operations_api` documentation for a complete list of operations. -Let's see how we can use this module and its operations. Import the module and -sub-module as +Import the module and sub-module as .. plot:: :format: doctest @@ -127,11 +127,10 @@ sub-module as >>> import mrsimulator.signal_processing as sp >>> import mrsimulator.signal_processing.apodization as apo -The signal processing is a series of operations 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. In the following -example, we show the application of a single operation—-a Gaussian line-broadening--using -the :class:`~mrsimulator.signal_processing.SignalProcessor` class, as follows, +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 @@ -146,10 +145,11 @@ the :class:`~mrsimulator.signal_processing.SignalProcessor` class, as follows, ... ) The required attribute of the ``SignalProcessor`` class, `operations`, holds the list of -operations that gets applied to the simulation dataset. In the above example, 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 +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:: @@ -162,9 +162,9 @@ transformation transforming the signal back to the frequency domain. Applying operation to the spectrum '''''''''''''''''''''''''''''''''' -To apply the above list of operations to the simulation data, use the +To apply the above list of operations to the simulation/input data, use the :meth:`~mrsimulator.signal_processing.SignalProcessor.apply_operations` method of the -instance as follows +``SignalProcessor`` instance as follows, .. plot:: :format: doctest @@ -174,8 +174,8 @@ instance as follows >>> 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 of the -signal processing operations. The plot of the processed signal is shown below. +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 @@ -234,7 +234,7 @@ Refer to the :ref:`config_simulator` section for further details. .. 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 are shown below +The plot of the spectra is shown below. .. plot:: :format: doctest @@ -247,7 +247,7 @@ The plot of the spectra are shown below >>> plt.show() # doctest: +SKIP Because the simulation is stored as a CSDM [#f1]_ object, each sub-spectrum is a -dependent-variable of the simulation dataset, sharing the same frequency dimension. +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 @@ -264,7 +264,7 @@ operations. ... operations=[ ... sp.IFFT(), # convert to time-domain ... apo.Gaussian(sigma=50, dep_var_indx=0), - ... apo.Exponential(Lambda=200, dep_var_indx=1), + ... apo.Exponential(FWHM=200, dep_var_indx=1), ... sp.FFT(), # convert to frequency-domain ... ] ... ) @@ -301,7 +301,7 @@ The plot of the processed spectrum is shown below. Serializing the operations list ------------------------------- -You may serialize the operations list using the +You may also serialize the operations list using the :meth:`~mrsimulator.signal_processing.SignalProcessor.to_dict_with_units` method, as follows @@ -315,14 +315,13 @@ method, as follows 'function': 'apodization', 'sigma': '50.0 Hz', 'type': 'Gaussian'}, - {'Lambda': '200.0 Hz', + {'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, 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/examples/1D_simulation/plot_0_Wollastonite.py b/examples/1D_simulation/plot_0_Wollastonite.py index f5ab6845f..21e02ce69 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 @@ -72,9 +74,16 @@ sim_wollastonite.run() # %% -# **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) + +# %% +# **Step 7** The plot of the simulation. 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..39d8879b7 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 @@ -59,9 +61,16 @@ sim_K2SO3.run() # %% -# **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) + +# %% +# **Step 7** The plot of the simulation. 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..0a1ce2c12 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 @@ -81,9 +83,16 @@ sim_coesite.run() # %% -# **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=145), sp.FFT()] +) +processed_data = post_sim.apply_operations(data=sim_coesite.methods[0].simulation) + +# %% +# **Step 7** The plot of the simulation. 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/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 58574c3ee..d67ff9a86 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -6,12 +6,12 @@ .. sectionauthor:: Maxwell C. Venetos """ # %% -# Often, after acquiring an NMR spectrum, we require some form of least-squares analysis -# to quantify our measurement. A typical recipe for any least-squares analysis comprises -# of two steps: +# 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 model parameters. +# - 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 @@ -116,14 +116,9 @@ import mrsimulator.signal_processing as sp import mrsimulator.signal_processing.apodization as apo -op_list = [ - sp.IFFT(), - apo.Exponential(Lambda=200), - sp.FFT(), - sp.Scale(factor=1), -] - -post_sim = sp.SignalProcessor(operations=op_list) +post_sim = sp.SignalProcessor( + operations=[sp.IFFT(), apo.Exponential(FWHM=200), sp.FFT(), sp.Scale(factor=0.3)] +) processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) # %% @@ -142,25 +137,24 @@ # # 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 LMFIT on its +# 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. , as follows +# 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=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="Lambda", value=post_sim.operations[1].Lambda) +params.add(name="FWHM", value=post_sim.operations[1].FWHM) params.add(name="factor", value=post_sim.operations[3].factor) # %% @@ -193,7 +187,7 @@ def minimization_function(params, sim, post_sim): # update the SignalProcessor parameter and apply line broadening. post_sim.operations[3].factor = values["factor"] - post_sim.operations[1].Lambda = values["Lambda"] + post_sim.operations[1].FWHM = values["FWHM"] processed_data = post_sim.apply_operations(sim.methods[0].simulation) # return the difference vector. @@ -202,14 +196,14 @@ def minimization_function(params, sim, post_sim): # %% # .. note:: -# To automate the fitting process, we provide a function to automatically parse +# 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 least-squares minimization -# '''''''''''''''''''''''''''''''''' +# 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* diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 95db62e22..247cd7dcc 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -144,7 +144,7 @@ factor = oxygen_experiment.max() op_list = [ sp.IFFT(), - apo.Exponential(Lambda=100), + apo.Exponential(FWHM=100), sp.FFT(), sp.Scale(factor=factor), ] diff --git a/src/mrsimulator/signal_processing/apodization.py b/src/mrsimulator/signal_processing/apodization.py index 65ea98a4d..20392b632 100644 --- a/src/mrsimulator/signal_processing/apodization.py +++ b/src/mrsimulator/signal_processing/apodization.py @@ -125,33 +125,76 @@ def operate(self, data): 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^{-\Lambda |{x}| \pi}, + f(x) = e^{-\Tau |x| \pi}, - where :math:`x` are the coordinates of the data dimension and :math:`\Lambda` is + 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 Lambda: The width parameter, :math:`\Lambda`. + 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) + >>> operation5 = apo.Exponential(FWHM=143.4, dim_indx=0, dep_var_indx=0) """ - Lambda: float = 0 + FWHM: float = 0 - property_unit_types: ClassVar = {"Lambda": ["time", "frequency"]} - property_default_units: ClassVar = {"Lambda": ["s", "Hz"]} - property_units: Dict = {"Lambda": "Hz"} + property_unit_types: ClassVar = {"FWHM": ["time", "frequency"]} + property_default_units: ClassVar = {"FWHM": ["s", "Hz"]} + property_units: Dict = {"FWHM": "Hz"} @staticmethod def fn(x, arg): @@ -165,6 +208,4 @@ def operate(self, data): dep_var: int. The index of the dependent variable to apply operation to """ - return self._operate( - data, fn=self.fn, prop_name="Lambda", prop_value=self.Lambda - ) + 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 index 39d1e79d2..db6792af3 100644 --- a/src/mrsimulator/signal_processing/tests/test_apodization.py +++ b/src/mrsimulator/signal_processing/tests/test_apodization.py @@ -28,7 +28,7 @@ PS_1 = [ sp.IFFT(dim_indx=0), - apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=0), + apo.Exponential(FWHM=200, dim_indx=0, dep_var_indx=0), sp.FFT(dim_indx=0), ] @@ -128,10 +128,10 @@ def test_scale_class(): def test_Exponential_class(): # direct initialization - a = apo.Exponential(Lambda=200, dim_indx=0, dep_var_indx=0) + a = apo.Exponential(FWHM=200, dim_indx=0, dep_var_indx=0) - assert a.Lambda == 200 - assert a.property_units == {"Lambda": "Hz"} + assert a.FWHM == 200 + assert a.property_units == {"FWHM": "Hz"} assert a.dim_indx == 0 assert a.dep_var_indx == 0 @@ -141,7 +141,7 @@ def test_Exponential_class(): assert dict_ == { "function": "apodization", "type": "Exponential", - "Lambda": "200.0 Hz", + "FWHM": "200.0 Hz", "dim_indx": 0, "dep_var_indx": 0, } diff --git a/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py index 3717b1ad4..915413210 100644 --- a/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py +++ b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py @@ -14,8 +14,8 @@ def test_01(): str2 = "spin_systems[0].sites[0].quadrupolar.Cq" str2_encoded = "ISO_0_SITES_0_quadrupolar_Cq" - assert sf._str_to_html(str1) == str1_encoded - assert sf._str_to_html(str2) == str2_encoded + assert sf._str_encode(str1) == str1_encoded + assert sf._str_encode(str2) == str2_encoded def test_02(): @@ -26,15 +26,15 @@ def test_02(): str2 = "spin_systems[0].sites[0].quadrupolar.Cq" str2_encoded = "ISO_0_SITES_0_quadrupolar_Cq" - assert sf._html_to_string(str1_encoded) == str1 - assert sf._html_to_string(str2_encoded) == str2 + 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(Lambda=100, dim_indx=0, dep_var_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), sp.FFT(dim_indx=0), sp.Scale(factor=10), ] @@ -52,7 +52,7 @@ def test_04(): # update_post_sim_from_LMFIT_params op_list = [ sp.IFFT(dim_indx=0), - apo.Exponential(Lambda=100, dim_indx=0, dep_var_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), sp.FFT(dim_indx=0), sp.Scale(factor=10), ] @@ -65,7 +65,7 @@ def test_04(): sf._update_post_sim_from_LMFIT_params(params, post_sim) - assert post_sim.operations[1].Lambda == 10 + assert post_sim.operations[1].FWHM == 10 assert post_sim.operations[3].factor == 1 @@ -84,7 +84,7 @@ def test_5(): op_list = [ sp.IFFT(dim_indx=0), - apo.Exponential(Lambda=100, dim_indx=0, dep_var_indx=0), + apo.Exponential(FWHM=100, dim_indx=0, dep_var_indx=0), sp.FFT(dim_indx=0), sp.Scale(factor=10), ] diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index a105ddbdd..fc08ab760 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -15,7 +15,7 @@ __email__ = "maxvenetos@gmail.com" START = "ISO_" -ENCRYPTION_PAIRS = [ +ENCODING_PAIRS = [ ["spin_systems[", START], ["].sites[", "_SITES_"], ["].isotropic_chemical_shift", "_isotropic_chemical_shift"], @@ -39,7 +39,7 @@ ] -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 @@ -51,12 +51,12 @@ def _str_to_html(my_string): Returns: String object. """ - for item in ENCRYPTION_PAIRS: + 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. @@ -67,7 +67,7 @@ def _html_to_string(my_string): Returns: String Object. """ - for item in ENCRYPTION_PAIRS: + for item in ENCODING_PAIRS: my_string = my_string.replace(*item[::-1]) return my_string @@ -106,13 +106,13 @@ def _traverse_dictionaries(dictionary, parent="spin_systems"): 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}]")) + name_list += _traverse_dictionaries(items, _str_encode(f"{parent}[{i}]")) return name_list @@ -138,7 +138,7 @@ def _post_sim_LMFIT_params(post_sim): temp_dict[f"{identifier}"] = arg elif isinstance(operation, apo.Exponential): identifier = f"operation_{i}_Exponential" - arg = operation.Lambda + arg = operation.FWHM temp_dict[f"{identifier}"] = arg elif isinstance(operation, sp.Scale): identifier = f"operation_{i}_Scale" @@ -160,7 +160,7 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): post_sim: SignalProcessor object """ temp_dict = {} - arg_dict = {"Gaussian": "sigma", "Exponential": "Lambda", "Scale": "factor"} + 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) @@ -228,13 +228,13 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): for items in temp_list: if "_eta" in items: params.add( - name=items, value=eval("sim." + _html_to_string(items)), min=0, max=1, + 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, @@ -242,12 +242,12 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): elif "abundance" in items: params.add( name=items, - value=eval("sim." + _html_to_string(items)) * abundance_scale, + value=eval("sim." + _str_decode(items)) * abundance_scale, min=0, max=100, ) else: - value = eval("sim." + _html_to_string(items)) + value = eval("sim." + _str_decode(items)) params.add(name=items, value=value) if post_sim is None: @@ -295,7 +295,7 @@ def LMFIT_min_function(params, sim, post_sim=None): values = params.valuesdict() for items in values: if "operation_" not in items: - nameString = "sim." + _html_to_string(items) + nameString = "sim." + _str_decode(items) executable = f"{nameString} = {values[items]}" exec(executable) elif "operation_" in items and post_sim is not None: From f514bdd82dc7d759de597b00d1a23c97ba8a9481 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sat, 25 Jul 2020 18:35:12 -0400 Subject: [PATCH 23/26] update csdmpy version to v0.3.1 update examples to include simulation and processed simulation. --- docs/requirements.rst | 2 +- docs/signal_processing.rst | 2 ++ environment-dev.yml | 2 +- environment.yml | 2 +- examples/1D_simulation/plot_0_Wollastonite.py | 12 ++++++++++-- examples/1D_simulation/plot_1_PotassiumSulfate.py | 13 +++++++++++-- examples/1D_simulation/plot_2_Coesite.py | 14 +++++++++++--- examples/Fitting/plot_2_mrsimFitExample_O17.py | 7 ++++++- requirements-dev.txt | 2 +- requirements.txt | 2 +- setup.py | 2 +- 11 files changed, 46 insertions(+), 14 deletions(-) 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 index 33dd04c76..80528146f 100644 --- a/docs/signal_processing.rst +++ b/docs/signal_processing.rst @@ -243,6 +243,7 @@ The plot of the spectra is shown below. >>> 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 @@ -294,6 +295,7 @@ The plot of the processed spectrum is shown below. >>> 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 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 21e02ce69..059ac1d58 100644 --- a/examples/1D_simulation/plot_0_Wollastonite.py +++ b/examples/1D_simulation/plot_0_Wollastonite.py @@ -22,6 +22,8 @@ # global plot configuration mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 2 + # %% # **Step 1** Create the sites. S29_1 = Site( @@ -73,6 +75,13 @@ # **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** Add post-simulation processing. post_sim = sp.SignalProcessor( @@ -80,8 +89,7 @@ ) processed_data = post_sim.apply_operations(data=sim_wollastonite.methods[0].simulation) -# %% -# **Step 7** The plot of the simulation. +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() diff --git a/examples/1D_simulation/plot_1_PotassiumSulfate.py b/examples/1D_simulation/plot_1_PotassiumSulfate.py index 39d8879b7..f3f62942b 100644 --- a/examples/1D_simulation/plot_1_PotassiumSulfate.py +++ b/examples/1D_simulation/plot_1_PotassiumSulfate.py @@ -21,6 +21,8 @@ # 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. S33 = Site( @@ -60,6 +62,14 @@ # **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** Add post-simulation processing. post_sim = sp.SignalProcessor( @@ -67,8 +77,7 @@ ) processed_data = post_sim.apply_operations(data=sim_K2SO3.methods[0].simulation) -# %% -# **Step 7** The plot of the simulation. +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() diff --git a/examples/1D_simulation/plot_2_Coesite.py b/examples/1D_simulation/plot_2_Coesite.py index 0a1ce2c12..c673d8a26 100644 --- a/examples/1D_simulation/plot_2_Coesite.py +++ b/examples/1D_simulation/plot_2_Coesite.py @@ -22,6 +22,8 @@ # global plot configuration mpl.rcParams["figure.figsize"] = [4.5, 3.0] +# sphinx_gallery_thumbnail_number = 2 + # %% # **Step 1** Create the sites. @@ -82,15 +84,21 @@ # **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** Add post-simulation processing. post_sim = sp.SignalProcessor( - operations=[sp.IFFT(), apo.Exponential(FWHM=30), apo.Gaussian(sigma=145), sp.FFT()] + 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) -# %% -# **Step 7** The plot of the simulation. +# The plot of the simulation after post-processing. ax = plt.subplot(projection="csdm") ax.plot(processed_data.real, color="black", linewidth=1) ax.invert_xaxis() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 247cd7dcc..ffc07a53f 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -57,6 +57,9 @@ # Convert the dimension coordinates from Hz to ppm. oxygen_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") +# Normalize the spectrum +oxygen_experiment /= oxygen_experiment.sum() + # plot of the dataset. ax = plt.subplot(projection="csdm") ax.plot(oxygen_experiment, color="black", linewidth=1) @@ -155,8 +158,10 @@ processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) ax = plt.subplot(projection="csdm") -ax.plot(processed_data.real, color="black", linewidth=1) +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.invert_xaxis() +plt.legend() plt.tight_layout() plt.show() 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", From 5e06de8aa33e6261905d79d40a76a745db09ce1c Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sat, 25 Jul 2020 20:10:52 -0400 Subject: [PATCH 24/26] examples update and format. --- examples/1D_simulation/plot_0_Wollastonite.py | 12 +++--- .../1D_simulation/plot_1_PotassiumSulfate.py | 12 +++--- examples/1D_simulation/plot_2_Coesite.py | 12 +++--- .../plot_1_mrsimFitExample_cuspidine.py | 15 ++++--- .../Fitting/plot_2_mrsimFitExample_O17.py | 42 +++++++++---------- 5 files changed, 49 insertions(+), 44 deletions(-) diff --git a/examples/1D_simulation/plot_0_Wollastonite.py b/examples/1D_simulation/plot_0_Wollastonite.py index 059ac1d58..2c7da0348 100644 --- a/examples/1D_simulation/plot_0_Wollastonite.py +++ b/examples/1D_simulation/plot_0_Wollastonite.py @@ -25,7 +25,7 @@ # 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 @@ -45,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 @@ -66,13 +66,13 @@ ) # %% -# **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. @@ -83,7 +83,7 @@ plt.show() # %% -# **Step 6** Add post-simulation processing. +# **Step 6:** Add post-simulation processing. post_sim = sp.SignalProcessor( operations=[sp.IFFT(), apo.Exponential(FWHM=70), sp.FFT()] ) diff --git a/examples/1D_simulation/plot_1_PotassiumSulfate.py b/examples/1D_simulation/plot_1_PotassiumSulfate.py index f3f62942b..2dd78d419 100644 --- a/examples/1D_simulation/plot_1_PotassiumSulfate.py +++ b/examples/1D_simulation/plot_1_PotassiumSulfate.py @@ -24,7 +24,7 @@ # 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", @@ -33,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 @@ -53,13 +53,13 @@ ) # %% -# **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. @@ -71,7 +71,7 @@ # %% -# **Step 6** Add post-simulation processing. +# **Step 6:** Add post-simulation processing. post_sim = sp.SignalProcessor( operations=[sp.IFFT(), apo.Exponential(FWHM=10), sp.FFT()] ) diff --git a/examples/1D_simulation/plot_2_Coesite.py b/examples/1D_simulation/plot_2_Coesite.py index c673d8a26..47c4c02f9 100644 --- a/examples/1D_simulation/plot_2_Coesite.py +++ b/examples/1D_simulation/plot_2_Coesite.py @@ -25,7 +25,7 @@ # 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( @@ -48,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 @@ -75,13 +75,13 @@ # 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. @@ -92,7 +92,7 @@ plt.show() # %% -# **Step 6** Add post-simulation processing. +# **Step 6:** Add post-simulation processing. post_sim = sp.SignalProcessor( operations=[sp.IFFT(), apo.Exponential(FWHM=30), apo.Gaussian(sigma=80), sp.FFT()] ) diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index d67ff9a86..3ba3c0faa 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -42,9 +42,13 @@ # convert the dimension coordinates from Hz to ppm synthetic_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") +# Normalize the spectrum +synthetic_experiment /= synthetic_experiment.max() + # 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() @@ -108,24 +112,25 @@ sim.methods[0].experiment = synthetic_experiment # %% -# **Step 5** simulate the spectrum. +# **Step 5:** simulate the spectrum. sim.run() # %% -# **Step 6** Create a SignalProcessor class and apply post simulation processing. +# **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=0.3)] + 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) # %% -# **Step 7** The plot the spectrum. We also plot the synthetic dataset for comparison. +# **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() @@ -224,10 +229,10 @@ def minimization_function(params, sim, post_sim): plt.plot(x, residual, alpha=0.5, label="Residual") 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() diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index ffc07a53f..df5aee93c 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -43,7 +43,6 @@ # Import the dataset # ------------------ # -# **Step 0** # Import the experimental data. In this example, we will import the data file serialized # with the CSDM file-format. We use the # `csdmpy `_ module to load the @@ -58,11 +57,12 @@ oxygen_experiment.dimensions[0].to("ppm", "nmr_frequency_ratio") # Normalize the spectrum -oxygen_experiment /= oxygen_experiment.sum() +oxygen_experiment /= oxygen_experiment.max() # 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() @@ -74,7 +74,7 @@ # Next, we will want to create a ``simulator`` object that we will use to fit to our # spectrum. We will start by creating a guess ``SpinSystem`` objects. # -# **Step 1** Create the guess sites. +# **Step 1:** Create the guess sites. O17_1 = Site( isotope="17O", isotropic_chemical_shift=60.0, # in ppm, @@ -88,11 +88,11 @@ ) # %% -# **Step 2** Create the spin systems for the guess sites. +# **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 fitting, you must create an +# **Step 3:** Create the method object. Note, when fitting, 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 @@ -125,41 +125,41 @@ method.experiment = oxygen_experiment # %% -# **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 = Simulator() sim.spin_systems = system_object sim.methods = [method] # %% -# **Step 5** Simulate the spectrum. +# **Step 5:** Simulate the spectrum. for iso in sim.spin_systems: # To avoid querying at every iteration we will save the relevant transition pathways iso.transition_pathways = method.get_transition_pathways(iso).tolist() sim.run() # %% -# **Step 6** Create the SignalProcessor class object and add to it the list of +# **Step 6:** Create the SignalProcessor class object and add to it the list of # operations that will be applied to the simulated spectrum. The generic list of # operations in NMR includes the line-broadening function. In the following, we add a # scaling and a Lorentzian line-broadening function. Here, the Lorentzian # line-broadening is defined as an exponential apodization operation sandwiched # between two Fourier transformations. -factor = oxygen_experiment.max() op_list = [ sp.IFFT(), - apo.Exponential(FWHM=100), + apo.Exponential(FWHM=40), sp.FFT(), - sp.Scale(factor=factor), + sp.Scale(factor=0.6), ] post_sim = sp.SignalProcessor(operations=op_list) # %% -# **Step 7** Process and plot the spectrum. +# **Step 7:** Process and plot the spectrum. processed_data = post_sim.apply_operations(data=sim.methods[0].simulation) 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() @@ -177,7 +177,7 @@ # attributes may vary. To simplify the parameter list creation we will use the # :func:`~mrsimulator.spectral_fitting.make_LMFIT_parameters` # -# **Step 8** Create a list of parameters to vary during fitting. +# **Step 8:** Create a list of parameters to vary during fitting. params = make_LMFIT_parameters(sim, post_sim) # %% @@ -189,25 +189,25 @@ params.pretty_print() # %% -# **Step 9** Perform the minimization. +# **Step 9:** Perform the minimization. minner = Minimizer(LMFIT_min_function, params, fcn_args=(sim, post_sim)) result = minner.minimize() report_fit(result) # %% -# **Step 10** Plot the fitted spectrum. +# **Step 10:** Plot the fitted spectrum. 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.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() From 80348bebad27f5cf223e76945d62b31732bca1e2 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 26 Jul 2020 03:02:25 -0400 Subject: [PATCH 25/26] update O17 least-squares example --- .../plot_1_mrsimFitExample_cuspidine.py | 2 +- .../Fitting/plot_2_mrsimFitExample_O17.py | 108 +++++++++--------- .../tests/test_spectral_fitting.py | 16 +-- src/mrsimulator/spectral_fitting.py | 24 ++-- 4 files changed, 75 insertions(+), 75 deletions(-) diff --git a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py index 3ba3c0faa..dae2855f3 100644 --- a/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py +++ b/examples/Fitting/plot_1_mrsimFitExample_cuspidine.py @@ -220,7 +220,7 @@ def minimization_function(params, sim, post_sim): print(fit_report(result)) # %% -# **The plot the fitted spectrum** +# The plot of the fit, measurement and the residuals is shown below. plt.figsize = (4, 3) x, y_data = synthetic_experiment.to_list() residual = result.residual diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index df5aee93c..00ef38b39 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -6,18 +6,16 @@ .. sectionauthor:: Maxwell C. Venetos """ # %% -# 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 relevant modules and presetting figure style and layout. import csdmpy as cp @@ -44,9 +42,8 @@ # ------------------ # # Import the experimental data. In this example, we will import the data file serialized -# with the CSDM file-format. We use the -# `csdmpy `_ module to load the -# dataset. +# with the CSDM file-format, using the +# `csdmpy `_ module. filename = "https://osu.box.com/shared/static/kfgt0jxgy93srsye9pofdnoha6qy58qf.csdf" oxygen_experiment = cp.load(filename) @@ -71,8 +68,8 @@ # Create a "fitting model" # ------------------------ # -# Next, we will want to create a ``simulator`` object that we will use to fit to our -# spectrum. We will start by creating a guess ``SpinSystem`` objects. +# 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( @@ -88,19 +85,19 @@ ) # %% -# **Step 2:** Create the spin systems for the guess sites. +# **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 fitting, 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 dimension metadata of the dimension. The remaining attribute values are set to -# the experimental conditions. +# **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 @@ -133,29 +130,22 @@ # %% # **Step 5:** Simulate the spectrum. for iso in sim.spin_systems: - # To avoid querying at every iteration we will save the relevant transition pathways + # To avoid querying at every iteration, save the relevant transition pathways info iso.transition_pathways = method.get_transition_pathways(iso).tolist() -sim.run() -# %% -# **Step 6:** Create the SignalProcessor class object and add to it the list of -# operations that will be applied to the simulated spectrum. The generic list of -# operations in NMR includes the line-broadening function. In the following, we add a -# scaling and a Lorentzian line-broadening function. Here, the Lorentzian -# line-broadening is defined as an exponential apodization operation sandwiched -# between two Fourier transformations. -op_list = [ - sp.IFFT(), - apo.Exponential(FWHM=40), - sp.FFT(), - sp.Scale(factor=0.6), -] -post_sim = sp.SignalProcessor(operations=op_list) +sim.run() # %% -# **Step 7:** Process and plot the spectrum. +# **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") @@ -169,33 +159,39 @@ # Least-squares minimization with LMFIT # ------------------------------------- # -# 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*. +# 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. # -# 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_LMFIT_parameters` -# -# **Step 8:** Create a list of parameters to vary during fitting. +# **Step 8:** Create a list of parameters. params = make_LMFIT_parameters(sim, post_sim) # %% +# 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 +# constraint on the fit by fixing the site abundances for the spin-systems at index # 1 and 2 to 50%. -params["ISO_0_abundance"].vary = False # fix the abundance +params["sys_0_abundance"].vary = False # fix the abundance params.pretty_print() # %% -# **Step 9:** Perform the minimization. +# **Step 9:** Perform least-squares minimization. We also provide a utility function, +# :func:`~mrsimulator.spectral_fitting.LMFIT_min_function`, for evaluating the +# difference vector between the simulation and experiment. You can 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() report_fit(result) # %% -# **Step 10:** Plot the fitted spectrum. +# **Step 10:** The plot of the fit, measurement and the residuals is shown below. plt.figsize = (4, 3) x, y_data = oxygen_experiment.to_list() residual = result.residual diff --git a/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py index 915413210..a9e1fb74c 100644 --- a/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py +++ b/src/mrsimulator/signal_processing/tests/test_spectral_fitting.py @@ -9,10 +9,10 @@ def test_01(): # str_to_html str1 = "spin_systems[0].sites[0].isotropic_chemical_shift" - str1_encoded = "ISO_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 = "ISO_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 @@ -21,10 +21,10 @@ def test_01(): def test_02(): # html_to_str str1 = "spin_systems[0].sites[0].isotropic_chemical_shift" - str1_encoded = "ISO_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 = "ISO_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 @@ -93,10 +93,10 @@ def test_5(): params = sf.make_LMFIT_parameters(sim, post_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, + "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, } diff --git a/src/mrsimulator/spectral_fitting.py b/src/mrsimulator/spectral_fitting.py index fc08ab760..4228a4b3d 100644 --- a/src/mrsimulator/spectral_fitting.py +++ b/src/mrsimulator/spectral_fitting.py @@ -14,10 +14,10 @@ __author__ = "Maxwell C Venetos" __email__ = "maxvenetos@gmail.com" -START = "ISO_" +START = "sys_" ENCODING_PAIRS = [ ["spin_systems[", START], - ["].sites[", "_SITES_"], + ["].sites[", "_site_"], ["].isotropic_chemical_shift", "_isotropic_chemical_shift"], ["].shielding_symmetric.", "_shielding_symmetric_"], ["].quadrupolar.", "_quadrupolar_"], @@ -194,8 +194,12 @@ def _update_post_sim_from_LMFIT_params(params, post_sim): def make_LMFIT_parameters(sim, post_sim=None, exclude_key=None): """ - Parses through the fitting parameter list to create LMFIT parameters used for - fitting. + 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. @@ -267,15 +271,15 @@ def make_LMFIT_parameters(sim, post_sim=None, exclude_key=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. From b11f63e828419f3f37cc02558197894443d94ccb Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Sun, 26 Jul 2020 14:11:24 -0400 Subject: [PATCH 26/26] update O17 fitting example. --- examples/Fitting/plot_2_mrsimFitExample_O17.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/Fitting/plot_2_mrsimFitExample_O17.py b/examples/Fitting/plot_2_mrsimFitExample_O17.py index 00ef38b39..f0d4b7d70 100644 --- a/examples/Fitting/plot_2_mrsimFitExample_O17.py +++ b/examples/Fitting/plot_2_mrsimFitExample_O17.py @@ -182,10 +182,11 @@ params.pretty_print() # %% -# **Step 9:** Perform least-squares minimization. We also provide a utility function, -# :func:`~mrsimulator.spectral_fitting.LMFIT_min_function`, for evaluating the -# difference vector between the simulation and experiment. You can use this function -# directly as the argument of the LMFIT Minimizer class, as follows, +# **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() report_fit(result)