From b3c53db5cb93b00063b1c3c967f36cf45721a614 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Wed, 11 Dec 2024 07:50:41 -0500 Subject: [PATCH 01/20] remove fista --- mrinversion/linear_model/_base_l1l2.py | 2 +- mrinversion/linear_model/fista_lasso.py | 6 ++- setup.py | 54 ++++++++++++------------- 3 files changed, 32 insertions(+), 30 deletions(-) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 58405f9..dca331c 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -534,7 +534,7 @@ def cv(l1, X, y, cv): y=y, scoring="neg_mean_squared_error", # 'neg_mean_absolute_error", cv=cv, - fit_params=fit_params, + params=fit_params, n_jobs=1, verbose=0, ) diff --git a/mrinversion/linear_model/fista_lasso.py b/mrinversion/linear_model/fista_lasso.py index aa3c383..06c9ff6 100644 --- a/mrinversion/linear_model/fista_lasso.py +++ b/mrinversion/linear_model/fista_lasso.py @@ -5,8 +5,10 @@ import numpy as np from mrinversion.linear_model._base_l1l2 import prepare_signal -from mrinversion.linear_model.fista import fista -from mrinversion.linear_model.fista import fista_cv + +# from mrinversion.linear_model.fista import fista +# from mrinversion.linear_model.fista import fista_cv +fista = fista_cv = None __author__ = "Deepansh Srivastava" CPU_COUNTS = os.cpu_count() diff --git a/setup.py b/setup.py index 50b21ac..54b2c6d 100644 --- a/setup.py +++ b/setup.py @@ -2,9 +2,10 @@ from os.path import dirname from os.path import join -from setuptools import find_packages -from numpy.distutils.core import setup -from numpy.distutils.core import Extension +from setuptools import find_packages, setup + +# from numpy.distutils.core import setup +# from numpy.distutils.core import Extension with open("mrinversion/__init__.py") as f: for line in f.readlines(): @@ -15,36 +16,36 @@ module_dir = dirname(abspath(__file__)) install_requires = [ - "numpy<2.0", + "numpy>2.0", "setuptools>=27.3", "csdmpy>=0.6", - "mrsimulator>=0.8.0rc0", + "mrsimulator>=1.0", "scikit-learn>=0.22", - "pydantic<2.0", + # "pydantic<2.0", ] setup_requires = ["setuptools>=27.3", "numpy"] extras = {"matplotlib": ["matplotlib>=3.0"]} -ext1 = Extension( - name="mrinversion.linear_model.fista.fista", - sources=["mrinversion/linear_model/fista/fista.f90"], - # f2py_options=['only:', 'subroutine_name', ':'], - extra_f90_compile_args=["-O3"], # '-fopenmp'], - libraries=["gomp", "blas"], - # f2py_options=['only:', 'fista', ':'], - language="f90", -) +# ext1 = Extension( +# name="mrinversion.linear_model.fista.fista", +# sources=["mrinversion/linear_model/fista/fista.f90"], +# # f2py_options=['only:', 'subroutine_name', ':'], +# extra_f90_compile_args=["-O3"], # '-fopenmp'], +# libraries=["gomp", "blas"], +# # f2py_options=['only:', 'fista', ':'], +# language="f90", +# ) -ext2 = Extension( - name="mrinversion.linear_model.fista.fista_cv", - sources=["mrinversion/linear_model/fista/fista_cv.f90"], - # f2py_options=['only:', 'subroutine_name', ':'], - extra_f90_compile_args=["-O3"], # '-fopenmp'], - libraries=["gomp", "blas"], - # f2py_options=['only:', 'fista', ':'], - language="f90", -) +# ext2 = Extension( +# name="mrinversion.linear_model.fista.fista_cv", +# sources=["mrinversion/linear_model/fista/fista_cv.f90"], +# # f2py_options=['only:', 'subroutine_name', ':'], +# extra_f90_compile_args=["-O3"], # '-fopenmp'], +# libraries=["gomp", "blas"], +# # f2py_options=['only:', 'fista', ':'], +# language="f90", +# ) setup( name="mrinversion", @@ -68,7 +69,7 @@ include_package_data=True, zip_safe=False, license="BSD-3-Clause", - ext_modules=[ext1, ext2], + # ext_modules=[ext1, ext2], classifiers=[ # Trove classifiers # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers @@ -78,11 +79,10 @@ "Development Status :: 4 - Beta", "License :: OSI Approved :: BSD License", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Topic :: Scientific/Engineering", ], ) From 65e19f8b00e0618197ce86562a7ea8216ac70de6 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 13 Dec 2024 15:52:08 -0500 Subject: [PATCH 02/20] add changes from PR 67 --- mrinversion/kernel/base.py | 11 +++- mrinversion/kernel/csa_aniso.py | 89 ++++++++++++++------------ mrinversion/kernel/nmr.py | 5 +- mrinversion/kernel/utils.py | 71 ++++++++++++++++++++ mrinversion/linear_model/_base_l1l2.py | 2 +- mrinversion/utils.py | 52 ++++++++++++++- 6 files changed, 184 insertions(+), 46 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 33a8ebb..3583bb8 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -4,6 +4,7 @@ import numpy as np from mrsimulator.method.lib import BlochDecaySpectrum +from .utils import _x_y_to_cq_eta_distribution from .utils import _x_y_to_zeta_eta_distribution __dimension_list__ = (cp.Dimension, cp.LinearDimension, cp.MonotonicDimension) @@ -162,6 +163,14 @@ def _get_zeta_eta(self, supersampling): ) return zeta, eta + def _get_cq_eta(self, supersampling): + """Return zeta and eta coordinates over x-y grid""" + + cq, eta = _x_y_to_cq_eta_distribution( + self.inverse_kernel_dimension, supersampling + ) + return cq, eta + def _check_csdm_dimension(dimensions, dimension_id): if not isinstance(dimensions, (list, *__dimension_list__)): @@ -201,4 +210,4 @@ def _check_dimension_type(dimensions, direction, dimension_quantity, kernel_type f"is required for the `{kernel_type}` kernel, instead got " f"`{item.quantity_name}` as the quantity name for the dimension at " f"index {i}." - ) + ) \ No newline at end of file diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index a27dcf8..14e8c7a 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -2,6 +2,7 @@ from mrsimulator import Simulator from mrsimulator import SpinSystem +from mrsimulator.method import Method from mrsimulator.method.lib import BlochDecaySpectrum from mrinversion.kernel.base import LineShape @@ -181,43 +182,51 @@ def __init__( ) -# class DAS(LineShape): -# def __init__( -# self, -# anisotropic_dimension, -# inverse_kernel_dimension, -# channel, -# magnetic_flux_density="9.4 T", -# rotor_angle="54.735 deg", -# rotor_frequency="600 Hz", -# number_of_sidebands=None, -# ): -# super().__init__( -# anisotropic_dimension, -# inverse_kernel_dimension, -# channel, -# magnetic_flux_density, -# rotor_angle, -# rotor_frequency, -# number_of_sidebands, -# # "DAS", -# ) - -# def kernel(self, supersampling): -# method = BlochDecayCentralTransitionSpectrum.parse_dict_with_units( -# self.method_args -# ) -# isotope = self.method_args["channels"][0] -# Cq, eta = self._get_zeta_eta(supersampling) -# spin_systems = [ -# SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) -# for cq_, e in zip(Cq, eta) -# ] - -# self.simulator.spin_systems = spin_systems -# self.simulator.methods = [method] -# self.simulator.run(pack_as_csdm=False) - -# amp = self.simulator.methods[0].simulation - -# return self._averaged_kernel(amp, supersampling) +class DAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # "DAS", + ) + + def kernel(self, supersampling): + # update method for DAS spectra events + das_event = dict( + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], + freq_contrib=["Quad2_4"], + ) + self.method_args["spectral_dimensions"][0]["events"] = [das_event] + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + Cq, eta = self._get_cq_eta(supersampling) + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + return self._averaged_kernel(amp, supersampling) \ No newline at end of file diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index d3bd1bf..451c1eb 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,7 +1,6 @@ +from mrinversion.kernel.csa_aniso import DAS # NOQA from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA -# from mrinversion.kernel.quad_aniso import MQMAS # NOQA - -# from mrinversion.kernel.lineshape import DAS # NOQA +# from mrinversion.kernel.quad_aniso import MQMAS # NOQA \ No newline at end of file diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index d60ee19..44603dc 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -129,6 +129,26 @@ def _x_y_to_zeta_eta_distribution(grid, supersampling): return _x_y_to_zeta_eta(x_mesh, y_mesh) +def _x_y_to_cq_eta_distribution(grid, supersampling): + """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" + x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) + y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) + + if x_coordinates.unit.physical_type == "frequency": + x_coordinates = x_coordinates.to("Hz").value + y_coordinates = y_coordinates.to("Hz").value + + elif x_coordinates.unit.physical_type == "dimensionless": + x_coordinates = x_coordinates.to("ppm").value + y_coordinates = y_coordinates.to("ppm").value + + x_mesh, y_mesh = np.meshgrid( + np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" + ) + + return _x_y_to_cq_eta(x_mesh, y_mesh) + + def _supersampled_coordinates(dimension, supersampling=1): r"""The coordinates along the dimension. @@ -171,3 +191,54 @@ def _supersampled_coordinates(dimension, supersampling=1): array[i::supersampling] = coordinates + (i - s2 + eo) * diff return array + + +def cq_eta_to_x_y(cq, eta): + r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the + following definition, + + .. math:: + \begin{array}{rl} + x &= C_q \sin\theta, \\ + y &= C_q \cos\theta + \end{array} {~~~~~~~~} + + where :math:`\theta = \frac{\pi}{2}\eta`. + + Args: + x: ndarray or list of floats. The coordinate x. + y: ndarray or list of floats. The coordinate y. + + Return: + A list of ndarrays. The first array holds the coordinate :math:`x`. The + second array holds the coordinates :math:`y`. + """ + cq = np.asarray(cq) + eta = np.asarray(eta) + + theta = np.pi * eta / 2.0 + x = np.zeros(cq.size) + y = np.zeros(cq.size) + + index = np.arange(len(cq)) + x[index] = cq[index] * np.cos(theta[index]) + y[index] = cq[index] * np.sin(theta[index]) + + return x.ravel(), y.ravel() + + +def _x_y_to_cq_eta(x, y): + """Same as def x_y_to_zeta_eta, but for ndarrays.""" + x = np.abs(x) + y = np.abs(y) + cq = np.sqrt(x**2 + y**2) # + offset + eta = np.ones(cq.shape) + # index = np.where(x > y) + index = np.arange(len(cq)) + # zeta[index] = -zeta[index] + eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) + + # index = np.where(x < y) + # eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index]) + + return cq.ravel(), eta.ravel() \ No newline at end of file diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index dca331c..81b7458 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -674,4 +674,4 @@ def prepare_signal(sig): scale = np.sqrt(np.mean(np.abs(s_) ** 2)) s_ = s_ / scale - return s_, scale + return s_, scale \ No newline at end of file diff --git a/mrinversion/utils.py b/mrinversion/utils.py index c4b55b9..bec9aeb 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -143,6 +143,56 @@ def get_polar_grids(ax, ticks=None, offset=0): ax.set_ylim(ylim) +def get_quadpolar_grids(ax, ticks=None, offset=0): + """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. + + Args: + ax: Matplotlib Axes. + ticks: Tick coordinates where radial grids are drawn. The value can be a list + or a numpy array. The default value is None. + offset: The grid is drawn at an offset away from the origin. + """ + ylim = ax.get_ylim() + xlim = ax.get_xlim() + if ticks is None: + x = np.asarray(ax.get_xticks()) + inc = x[1] - x[0] + size = x.size + x = np.arange(size + 5) * inc + else: + x = np.asarray(ticks) + + lw = 0.3 + # t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) + # t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) + + # ax.add_artist(t1) + # ax.add_artist(t2) + for x_ in x: + if x_ - offset > 0: + ax.add_artist( + plt.Circle( + (0, 0), + x_ - offset, + fill=False, + color="k", + linestyle="--", + linewidth=lw, + alpha=0.5, + ) + ) + + angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 2.0) + # angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) + for ang_ in angle1: + ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) + # for ang_ in angle2: + # ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) + # ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + def plot_3d( ax, csdm_objects, @@ -448,4 +498,4 @@ def plot_3d_( ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) + ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) \ No newline at end of file From 4e1a7a81f29b28d48f11f34de64b4233c1d13fc1 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Thu, 30 Jan 2025 15:59:00 -0500 Subject: [PATCH 03/20] try to add masked array option --- mrinversion/kernel/base.py | 11 +- mrinversion/kernel/csa_aniso.py | 5 +- mrinversion/linear_model/linear_inversion.py | 5 + mrinversion/utils.py | 120 +++++++++++++++++++ 4 files changed, 137 insertions(+), 4 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 3583bb8..4e8b3b5 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -2,6 +2,7 @@ import csdmpy as cp import numpy as np +import numpy.ma as ma from mrsimulator.method.lib import BlochDecaySpectrum from .utils import _x_y_to_cq_eta_distribution @@ -52,7 +53,7 @@ def __init__(self, kernel_dimension, inverse_kernel_dimension, n_dir, n_inv): self.inverse_kernel_dimension = inverse_kernel_dimension self.inverse_dimension = self.inverse_kernel_dimension - def _averaged_kernel(self, amp, supersampling, xy_grid=True): + def _averaged_kernel(self, amp, supersampling, xy_grid=True, mask_kernel=False): """Return the kernel by averaging over the supersampled grid cells.""" shape = () inverse_kernel_dimension = self.inverse_kernel_dimension @@ -79,7 +80,13 @@ def _averaged_kernel(self, amp, supersampling, xy_grid=True): section_ = deepcopy(section) section_[i] = 0 K[tuple(section_)] /= 2.0 - + print(f'K-shape: {K.shape}') + if mask_kernel: + mask = np.zeros(K.shape) + for i in range(len(K)): + mask[i][len(K[0])-i:] += 1 + K = ma.masked_array(K, mask=mask) + inv_size = np.asarray([item.count for item in inverse_kernel_dimension]).prod() K = K.reshape(inv_size, self.kernel_dimension.count).T return K diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index 14e8c7a..49f752f 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -204,7 +204,7 @@ def __init__( # "DAS", ) - def kernel(self, supersampling): + def kernel(self, supersampling, mask_kernel=False): # update method for DAS spectra events das_event = dict( transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], @@ -215,6 +215,7 @@ def kernel(self, supersampling): method = Method.parse_dict_with_units(self.method_args) isotope = self.method_args["channels"][0] Cq, eta = self._get_cq_eta(supersampling) + Cq, eta = self._get_zeta_eta(supersampling) spin_systems = [ SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) for cq_, e in zip(Cq, eta) @@ -229,4 +230,4 @@ def kernel(self, supersampling): amp = sim.methods[0].simulation.real - return self._averaged_kernel(amp, supersampling) \ No newline at end of file + return self._averaged_kernel(amp, supersampling, mask_kernel=mask_kernel) \ No newline at end of file diff --git a/mrinversion/linear_model/linear_inversion.py b/mrinversion/linear_model/linear_inversion.py index 363644f..9900031 100644 --- a/mrinversion/linear_model/linear_inversion.py +++ b/mrinversion/linear_model/linear_inversion.py @@ -1,4 +1,5 @@ import numpy as np +# import numpy.ma as ma __author__ = "Deepansh J. Srivastava" __email__ = "srivastava.89@osu.edu" @@ -24,6 +25,10 @@ def find_optimum_singular_value(s): def TSVD(K): + # print(f'K: {K}') + # if ma.is_masked(K): + # U, S, VT = np.linalg.svd(K.filled(0), full_matrices=False) + # else: U, S, VT = np.linalg.svd(K, full_matrices=False) r = find_optimum_singular_value(S) return U, S, VT, r diff --git a/mrinversion/utils.py b/mrinversion/utils.py index bec9aeb..33d7c32 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -93,6 +93,126 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): return csdm_new +def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): + """Convert the three-dimensional p(iso, Cq, eta) to p(iso, xq, yq) tensor + distribution. + + Args + ---- + + csdm_object: CSDM + A CSDM object containing the 3D p(iso, Cq, eta) distribution. + xq: CSDM.Dimension + A CSDM dimension object describing the xq dimension. + yq: CSDM.Dimension + A CSDM dimension object describing the yq dimension. + n: int + An integer used in linear interpolation of the data. The default is 5. + """ + [ + item.to("ppm", "nmr_frequency_ratio") + for item in csdm_object.x + if item.origin_offset != 0 + ] + data = csdm_object.y[0].components[0] + print(data) + # csdm_object.plot() + extra_dims = 1 + if len(csdm_object.x) > 2: + extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) + data.shape = (extra_dims, data.shape[-2], data.shape[-1]) + + reg_cq, reg_eta = (csdm_object.x[i].coordinates.value for i in range(2)) + d_cq = reg_cq[1] - reg_cq[0] + d_eta = reg_eta[1] - reg_eta[0] + sol = np.zeros((extra_dims, xq.count, yq.count)) + + bins = [xq.count, yq.count] + dx = xq.increment.value / 2 + dy = yq.increment.value / 2 + range_ = [ + [xq.coordinates[0].value - dx, xq.coordinates[-1].value + dx], + [yq.coordinates[0].value - dy, yq.coordinates[-1].value + dy], + ] + # print(range_) + avg_range_cq = (np.arange(n) - (n - 1) / 2) * d_cq / n + avg_range_eta = (np.arange(n) - (n - 1) / 2) * d_eta / n + # print(avg_range_cq) + # print(avg_range_eta) + for i,cq_item in enumerate(avg_range_cq): + for j,eta_item in enumerate(avg_range_eta): + # print(i,j) + cq__ = (reg_cq + cq_item) + # print(np.abs(cq__).shape) + eta__ = np.abs(reg_eta + eta_item) + cq_, eta_ = np.meshgrid(cq__, eta__) + # print(cq_) + # print() + cq_ = cq_.ravel() + eta_ = eta_.ravel() + # print(cq_.shape) + # print(eta_.shape) + # print(eta_) + # print() + + theta = np.ones(eta_.shape)*1e10 + # print(theta.shape) + + index = np.where(cq_ > 0) + # print(f'pos index: {index}') + theta[index] = np.pi/2 * (1-eta_[index] / 2.0) + + index = np.where(cq_ <= 0) + # print(f'neg index: {index}') + + theta[index] = np.pi/4.0 * eta_[index] + # print(cq_) + # print(np.where(theta * 180/np.pi==1e10)) + # print(list(zip(cq_, theta * 180/np.pi))) + + + xq_grid = np.abs(cq_) * np.sin(theta) + yq_grid = np.abs(cq_) * np.cos(theta) + + # index = np.arange(xq.count) + # print(f'index: {index}') + # print(data[0].ravel()) + # print(data[0].shape) + ################eta_grid = (2 / np.pi) * np.arctan(y_ / x_) + # print(f'extra dims: {extra_dims}') + for i in range(extra_dims): + weight = deepcopy(data[i]).ravel() + # weight[index] /= 2 + # np.append(weight, weight[index]) + # plt.plot(weight) + # plt.show() + # print(np.where(weight>0.1)) + # print(range_) + sol_, _, _ = np.histogram2d( + xq_grid, yq_grid, weights=weight, bins=bins, range=range_ + ) + sol[i] += sol_ + # print(sol) + # print(np.where(sol>0.1)) + # plt.imshow(sol[0,:,:]) + # plt.imshow(sol_) + # print(xq.coordinates.shape) + # print(yq.coordinates.shape) + # plt.plot(weight) + # plt.show() + + sol /= n * n + # plt.plot(weight) + # plt.show() + del xq_grid, yq_grid, index, cq_, eta_, avg_range_cq, avg_range_eta + csdm_new = cp.as_csdm(np.squeeze(sol)) + csdm_new.x[0] = yq + csdm_new.x[1] = xq + if len(csdm_object.x) > 2: + csdm_new.x[2] = csdm_object.x[2] + return csdm_new + + def get_polar_grids(ax, ticks=None, offset=0): """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. From 1b021a06347057acafa17fb714f5fdf008bad7a8 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 9 May 2025 22:52:08 -0400 Subject: [PATCH 04/20] add functions --- mrinversion/kernel/base.py | 33 ++++++---- mrinversion/kernel/csa_aniso.py | 35 ++++++++--- mrinversion/kernel/utils.py | 103 ++++++++++++++++++++++++++++---- mrinversion/utils.py | 11 ++-- 4 files changed, 147 insertions(+), 35 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 4e8b3b5..9c48b2b 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -63,15 +63,18 @@ def _averaged_kernel(self, amp, supersampling, xy_grid=True, mask_kernel=False): for item in inverse_kernel_dimension[::-1]: shape += (item.count, supersampling) shape += (self.kernel_dimension.count,) - + # print(f'shape: {shape}') K = amp.reshape(shape) inv_len = len(inverse_kernel_dimension) axes = tuple(2 * i + 1 for i in range(inv_len)) K = K.mean(axis=axes) + # print(f'K shape: {K.shape}') if xy_grid: section = [*[0 for _ in range(inv_len)], slice(None, None, None)] + # print(f'section: {section}') + # print(f'K[tuple(section)]: {K[tuple(section)]}') K /= K[tuple(section)].sum() section = [slice(None, None, None) for _ in range(inv_len + 1)] @@ -80,12 +83,12 @@ def _averaged_kernel(self, amp, supersampling, xy_grid=True, mask_kernel=False): section_ = deepcopy(section) section_[i] = 0 K[tuple(section_)] /= 2.0 - print(f'K-shape: {K.shape}') - if mask_kernel: - mask = np.zeros(K.shape) - for i in range(len(K)): - mask[i][len(K[0])-i:] += 1 - K = ma.masked_array(K, mask=mask) + # print(f'K-shape: {K.shape}') + # if mask_kernel: + # mask = np.zeros(K.shape) + # for i in range(len(K)): + # mask[i][len(K[0])-i:] += 1 + # K = ma.masked_array(K, mask=mask) inv_size = np.asarray([item.count for item in inverse_kernel_dimension]).prod() K = K.reshape(inv_size, self.kernel_dimension.count).T @@ -162,13 +165,19 @@ def __init__( if number_of_sidebands is None: self.number_of_sidebands = dim.count - def _get_zeta_eta(self, supersampling): + def _get_zeta_eta(self, supersampling, eta_bound=1): """Return zeta and eta coordinates over x-y grid""" - zeta, eta = _x_y_to_zeta_eta_distribution( - self.inverse_kernel_dimension, supersampling - ) - return zeta, eta + if eta_bound == 1: + zeta, eta = _x_y_to_zeta_eta_distribution( + self.inverse_kernel_dimension, supersampling, eta_bound + ) + return zeta, eta + else: + zeta, eta, abundances= _x_y_to_zeta_eta_distribution( + self.inverse_kernel_dimension, supersampling, eta_bound + ) + return zeta, eta, abundances def _get_cq_eta(self, supersampling): """Return zeta and eta coordinates over x-y grid""" diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index 49f752f..ef0c78d 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -204,7 +204,7 @@ def __init__( # "DAS", ) - def kernel(self, supersampling, mask_kernel=False): + def kernel(self, supersampling, mask_kernel=False, return_as_obj=False, eta_bound = 1): # update method for DAS spectra events das_event = dict( transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], @@ -214,12 +214,28 @@ def kernel(self, supersampling, mask_kernel=False): method = Method.parse_dict_with_units(self.method_args) isotope = self.method_args["channels"][0] - Cq, eta = self._get_cq_eta(supersampling) - Cq, eta = self._get_zeta_eta(supersampling) - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] + # Cq, eta = self._get_cq_eta(supersampling) + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + # Cq = list(Cq) + # eta = list(eta) + # print(f'len Cq: {len(Cq)}; len eta: {len(eta)}') + # for i, (_, eta_) in enumerate(zip(Cq, eta)): + # if eta_ > eta_bound: + # del Cq[i], eta[i] + # print(f'len Cq: {len(Cq)}; len eta: {len(eta)}') + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] sim = Simulator() sim.config.number_of_sidebands = self.number_of_sidebands sim.config.decompose_spectrum = "spin_system" @@ -229,5 +245,8 @@ def kernel(self, supersampling, mask_kernel=False): sim.run(pack_as_csdm=False) amp = sim.methods[0].simulation.real - + # print(Cq) + # print(eta) + if return_as_obj: + kernel = self._averaged_kernel(amp, supersampling, mask_kernel=mask_kernel) return self._averaged_kernel(amp, supersampling, mask_kernel=mask_kernel) \ No newline at end of file diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index 44603dc..e19fe32 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -49,20 +49,101 @@ def x_y_to_zeta_eta(x, y): return zeta * x_unit, eta -def _x_y_to_zeta_eta(x, y): +def _x_y_to_zeta_eta(x, y, eta_bound=1): """Same as def x_y_to_zeta_eta, but for ndarrays.""" - x = np.abs(x) - y = np.abs(y) + x = np.abs(x).ravel() + y = np.abs(y).ravel() + # print(x/1e6) + # print(y/1e6) + n_lost = 0 zeta = np.sqrt(x**2 + y**2) # + offset + abundances = np.ones(zeta.shape) eta = np.ones(zeta.shape) - index = np.where(x > y) - zeta[index] = -zeta[index] - eta[index] = (4.0 / np.pi) * np.arctan(y[index] / x[index]) + eta = eta.tolist() + # print(zeta/1e6) + zeta = zeta.tolist() + # x.tolist() + # y.tolist() + + # print(f'eta length pre: {len(eta)}') + # print(eta.shape) + del_these = [] + for index,_ in enumerate(x): + if x[index] > y[index]: + zeta[index] = - zeta[index] + this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) + if this_eta < eta_bound: + eta[index] = this_eta + else: + # del_these.append(index) + # del zeta[index], eta[index] + # n_lost += 1 + abundances[index] = 0 + elif x[index] < y[index]: + this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) + if this_eta < eta_bound: + eta[index] = this_eta + else: + # del zeta[index], eta[index] + # del_these.append(index) + # n_lost += 1 + abundances[index] = 0 + elif x[index] == y[index] and eta_bound < 1: + # del_these.append(index) + # n_lost += 1 + abundances[index] = 0 + + # print(np.asarray(eta)) + # print(del_these) + # print(f'eta length post: {len(eta)}') + # print(f'n lost: {n_lost}') + if eta_bound ==1: + return np.asarray(zeta), np.asarray(eta) + else: + for idx in del_these[::-1]: + del zeta[idx], eta[idx] + # print(np.asarray(eta)) + return np.asarray(zeta), np.asarray(eta), abundances + + + + + +# def _x_y_to_zeta_eta(x, y, eta_bound=1): +# """Same as def x_y_to_zeta_eta, but for ndarrays.""" +# x = np.abs(x) +# y = np.abs(y) +# # print(x) +# # print(y) +# n_lost = 0 +# zeta = np.sqrt(x**2 + y**2) # + offset +# # print(zeta) +# eta = np.ones(zeta.shape)*100 +# # print(eta.shape) +# for index,_ in np.ndenumerate(x): +# if x[index] > y[index]: +# zeta[index] = - zeta[index] +# this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) +# if this_eta < eta_bound: +# eta[index] = this_eta +# else: +# n_lost += 1 +# elif x[index] < y[index]: +# this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) +# if this_eta < eta_bound: +# eta[index] = this_eta +# else: +# n_lost += 1 +# print(eta) + +# if eta_bound ==1: +# return zeta.ravel(), eta.ravel() +# else: +# return zeta.ravel(), eta.ravel(), n_lost + + - index = np.where(x < y) - eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index]) - return zeta.ravel(), eta.ravel() def zeta_eta_to_x_y(zeta, eta): @@ -109,7 +190,7 @@ def zeta_eta_to_x_y(zeta, eta): return x.ravel(), y.ravel() -def _x_y_to_zeta_eta_distribution(grid, supersampling): +def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound): """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) @@ -126,7 +207,7 @@ def _x_y_to_zeta_eta_distribution(grid, supersampling): np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" ) - return _x_y_to_zeta_eta(x_mesh, y_mesh) + return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound) def _x_y_to_cq_eta_distribution(grid, supersampling): diff --git a/mrinversion/utils.py b/mrinversion/utils.py index 33d7c32..c8ce5d6 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -31,12 +31,12 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): if item.origin_offset != 0 ] data = csdm_object.y[0].components[0] - + # print(f'data max: {data.max()}') extra_dims = 1 if len(csdm_object.x) > 2: extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) data.shape = (extra_dims, data.shape[-2], data.shape[-1]) - + reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) dx = reg_x[1] - reg_x[0] dy = reg_y[1] - reg_y[0] @@ -47,7 +47,7 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): d_eta = eta.increment.value / 2 range_ = [ [zeta.coordinates[0].value - d_zeta, zeta.coordinates[-1].value + d_zeta], - [eta.coordinates[0] - d_eta, eta.coordinates[-1] + d_eta], + [eta.coordinates[0].value - d_eta, eta.coordinates[-1].value + d_eta], ] avg_range_x = (np.arange(n) - (n - 1) / 2) * dx / n @@ -75,13 +75,16 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): np.append(eta, np.ones(index[0].size)) for i in range(extra_dims): weight = deepcopy(data[i]).ravel() + # print(f'weight max: {weight.max()}') + # print(f'range: {range_}') weight[index] /= 2 np.append(weight, weight[index]) sol_, _, _ = np.histogram2d( zeta_grid, eta_grid, weights=weight, bins=bins, range=range_ ) sol[i] += sol_ - + # print(f'sol max: {sol.max()}') + # print() sol /= n * n del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y From f72eb13972119584923e3fc3f3cb4718678aa0cd Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava <21365911+deepanshs@users.noreply.github.com> Date: Mon, 2 Jun 2025 05:45:53 -0400 Subject: [PATCH 05/20] remove commented fista import --- mrinversion/linear_model/fista_lasso.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/mrinversion/linear_model/fista_lasso.py b/mrinversion/linear_model/fista_lasso.py index 5c274b0..d6a78ee 100644 --- a/mrinversion/linear_model/fista_lasso.py +++ b/mrinversion/linear_model/fista_lasso.py @@ -5,10 +5,8 @@ import numpy as np from mrinversion.linear_model._base_l1l2 import prepare_signal - -# from mrinversion.linear_model.fista import fista -# from mrinversion.linear_model.fista import fista_cv -fista = fista_cv = None +from mrinversion.linear_model.fista import fista +from mrinversion.linear_model.fista import fista_cv __author__ = "Deepansh Srivastava" CPU_COUNTS = os.cpu_count() From cc18d42fe08b62cf4453904d58ea04b403171f43 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Mon, 23 Jun 2025 23:27:24 -0400 Subject: [PATCH 06/20] static limit mqmas kernel --- mrinversion/kernel/csa_aniso.py | 255 ++++++++++++++++++++++++++++++-- mrinversion/kernel/nmr.py | 3 + 2 files changed, 243 insertions(+), 15 deletions(-) diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index ef0c78d..5a391da 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -1,8 +1,12 @@ from copy import deepcopy +import csdmpy as cp +import numpy as np +from mrsimulator.utils import get_spectral_dimensions + from mrsimulator import Simulator from mrsimulator import SpinSystem -from mrsimulator.method import Method +from mrsimulator.method import Method, SpectralEvent from mrsimulator.method.lib import BlochDecaySpectrum from mrinversion.kernel.base import LineShape @@ -204,7 +208,7 @@ def __init__( # "DAS", ) - def kernel(self, supersampling, mask_kernel=False, return_as_obj=False, eta_bound = 1): + def kernel(self, supersampling, eta_bound = 1): # update method for DAS spectra events das_event = dict( transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], @@ -214,18 +218,80 @@ def kernel(self, supersampling, mask_kernel=False, return_as_obj=False, eta_boun method = Method.parse_dict_with_units(self.method_args) isotope = self.method_args["channels"][0] - # Cq, eta = self._get_cq_eta(supersampling) + + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + return self._averaged_kernel(amp, supersampling) + + +class MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # "DAS", + ) + + def kernel(self, supersampling, eta_bound = 1): + # update method for DAS spectra events + mqmas_events = [ + dict( + fraction= -9/50, + transition_queries=[{"ch1": {"P": [-3], "D": [0]}}] + ), + dict( + fraction= 27/50, + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}] + ), + ] + + self.method_args["spectral_dimensions"][0]["events"] = mqmas_events + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + if eta_bound == 1: Cq, eta = self._get_zeta_eta(supersampling, eta_bound) else: Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - # Cq = list(Cq) - # eta = list(eta) - # print(f'len Cq: {len(Cq)}; len eta: {len(eta)}') - # for i, (_, eta_) in enumerate(zip(Cq, eta)): - # if eta_ > eta_bound: - # del Cq[i], eta[i] - # print(f'len Cq: {len(Cq)}; len eta: {len(eta)}') + if eta_bound == 1: spin_systems = [ SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) @@ -245,8 +311,167 @@ def kernel(self, supersampling, mask_kernel=False, return_as_obj=False, eta_boun sim.run(pack_as_csdm=False) amp = sim.methods[0].simulation.real - # print(Cq) - # print(eta) - if return_as_obj: - kernel = self._averaged_kernel(amp, supersampling, mask_kernel=mask_kernel) - return self._averaged_kernel(amp, supersampling, mask_kernel=mask_kernel) \ No newline at end of file + print(f'amp shape: {amp.shape}') + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMASnodist(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # + # "DAS", + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound = 1): + import sys + sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') + # import src.processing as smproc + import src.simulation as smsim + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + # print(f'Cq: {Cq}') + # print(f'eta: {eta}') + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + # print(self.anisotropic_dimension) + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + # print(obj.x[0]) + # print(self.anisotropic_dimension) + # print(spec_dim) + + amp = np.asarray([smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type='c0_c4', + contribs='c0_c4', + return_array=True, + distorted=False) for mysys in spin_systems]) + # sim = Simulator() + # sim.config.number_of_sidebands = self.number_of_sidebands + # sim.config.decompose_spectrum = "spin_system" + + # sim.spin_systems = spin_systems + # sim.methods = [method] + # sim.run(pack_as_csdm=False) + # obj_pre = cp.CSDM( + # dimensions=[ + # aniso + # ] + # ) + # amp = sim.methods[0].simulation.real + # print(f'amp shape: {amp.shape}') + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # + # "DAS", + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound = 1): + import sys + sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') + # import src.processing as smproc + import src.simulation as smsim + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + print(f'Cq: {Cq}') + print(f'eta: {eta}') + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + # print(self.anisotropic_dimension) + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray([smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type='c0_c4', + contribs='c0_c4', + return_array=True, + distorted=True) for mysys in spin_systems]) + # sim = Simulator() + # sim.config.number_of_sidebands = self.number_of_sidebands + # sim.config.decompose_spectrum = "spin_system" + + # sim.spin_systems = spin_systems + # sim.methods = [method] + # sim.run(pack_as_csdm=False) + + # amp = sim.methods[0].simulation.real + print(f'amp shape: {amp.shape}') + return self._averaged_kernel(amp, supersampling) + \ No newline at end of file diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index 451c1eb..32d92d3 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,4 +1,7 @@ from mrinversion.kernel.csa_aniso import DAS # NOQA +from mrinversion.kernel.csa_aniso import MQMAS # NOQA +from mrinversion.kernel.csa_aniso import SL_MQMAS # NOQA +from mrinversion.kernel.csa_aniso import SL_MQMASnodist from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA From 9a8ea8bddc22aa7989afce6efe700d4d1370e0d1 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:22:46 -0400 Subject: [PATCH 07/20] add CV minimization and xygrid parameter to SmoothLasso and SmoothLassoCV --- mrinversion/kernel/base.py | 7 +- mrinversion/kernel/csa_aniso.py | 289 ----------------------- mrinversion/kernel/nmr.py | 8 +- mrinversion/kernel/quad_aniso.py | 266 +++++++++++++++++++++ mrinversion/kernel/utils.py | 186 +++++---------- mrinversion/linear_model/_base_l1l2.py | 139 +++++++++-- mrinversion/linear_model/smooth_lasso.py | 24 ++ mrinversion/linear_model/utils.py | 199 ++++++++++++++++ 8 files changed, 680 insertions(+), 438 deletions(-) create mode 100644 mrinversion/linear_model/utils.py diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 9c48b2b..00eb4ea 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -165,17 +165,16 @@ def __init__( if number_of_sidebands is None: self.number_of_sidebands = dim.count - def _get_zeta_eta(self, supersampling, eta_bound=1): + def _get_zeta_eta(self, supersampling, eta_bound=1,calc_pos=False): """Return zeta and eta coordinates over x-y grid""" - - if eta_bound == 1: + if eta_bound == 1 and not calc_pos: zeta, eta = _x_y_to_zeta_eta_distribution( self.inverse_kernel_dimension, supersampling, eta_bound ) return zeta, eta else: zeta, eta, abundances= _x_y_to_zeta_eta_distribution( - self.inverse_kernel_dimension, supersampling, eta_bound + self.inverse_kernel_dimension, supersampling, eta_bound, calc_pos ) return zeta, eta, abundances diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index 5a391da..609aee5 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -186,292 +186,3 @@ def __init__( ) -class DAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # "DAS", - ) - - def kernel(self, supersampling, eta_bound = 1): - # update method for DAS spectra events - das_event = dict( - transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], - freq_contrib=["Quad2_4"], - ) - self.method_args["spectral_dimensions"][0]["events"] = [das_event] - - method = Method.parse_dict_with_units(self.method_args) - isotope = self.method_args["channels"][0] - - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - sim = Simulator() - sim.config.number_of_sidebands = self.number_of_sidebands - sim.config.decompose_spectrum = "spin_system" - - sim.spin_systems = spin_systems - sim.methods = [method] - sim.run(pack_as_csdm=False) - - amp = sim.methods[0].simulation.real - - return self._averaged_kernel(amp, supersampling) - - -class MQMAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # "DAS", - ) - - def kernel(self, supersampling, eta_bound = 1): - # update method for DAS spectra events - mqmas_events = [ - dict( - fraction= -9/50, - transition_queries=[{"ch1": {"P": [-3], "D": [0]}}] - ), - dict( - fraction= 27/50, - transition_queries=[{"ch1": {"P": [-1], "D": [0]}}] - ), - ] - - self.method_args["spectral_dimensions"][0]["events"] = mqmas_events - - method = Method.parse_dict_with_units(self.method_args) - isotope = self.method_args["channels"][0] - - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - sim = Simulator() - sim.config.number_of_sidebands = self.number_of_sidebands - sim.config.decompose_spectrum = "spin_system" - - sim.spin_systems = spin_systems - sim.methods = [method] - sim.run(pack_as_csdm=False) - - amp = sim.methods[0].simulation.real - print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - - -class SL_MQMASnodist(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - exp_dict, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # - # "DAS", - ) - self.exp_dict = exp_dict - self.anisotropic_dimension = anisotropic_dimension - - def kernel(self, supersampling, eta_bound = 1): - import sys - sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') - # import src.processing as smproc - import src.simulation as smsim - # import src.fitting as smfit - - isotope = self.method_args["channels"][0] - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - # print(f'Cq: {Cq}') - # print(f'eta: {eta}') - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - # print(self.anisotropic_dimension) - - obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) - spec_dim = get_spectral_dimensions(obj) - # print(obj.x[0]) - # print(self.anisotropic_dimension) - # print(spec_dim) - - amp = np.asarray([smsim.simulate_onesite_lineshape( - self.exp_dict, - mysys, - spec_dim[0], - input_type='c0_c4', - contribs='c0_c4', - return_array=True, - distorted=False) for mysys in spin_systems]) - # sim = Simulator() - # sim.config.number_of_sidebands = self.number_of_sidebands - # sim.config.decompose_spectrum = "spin_system" - - # sim.spin_systems = spin_systems - # sim.methods = [method] - # sim.run(pack_as_csdm=False) - # obj_pre = cp.CSDM( - # dimensions=[ - # aniso - # ] - # ) - # amp = sim.methods[0].simulation.real - # print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - - -class SL_MQMAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - exp_dict, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # - # "DAS", - ) - self.exp_dict = exp_dict - self.anisotropic_dimension = anisotropic_dimension - - def kernel(self, supersampling, eta_bound = 1): - import sys - sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') - # import src.processing as smproc - import src.simulation as smsim - # import src.fitting as smfit - - isotope = self.method_args["channels"][0] - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - print(f'Cq: {Cq}') - print(f'eta: {eta}') - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - # print(self.anisotropic_dimension) - - obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) - spec_dim = get_spectral_dimensions(obj) - - amp = np.asarray([smsim.simulate_onesite_lineshape( - self.exp_dict, - mysys, - spec_dim[0], - input_type='c0_c4', - contribs='c0_c4', - return_array=True, - distorted=True) for mysys in spin_systems]) - # sim = Simulator() - # sim.config.number_of_sidebands = self.number_of_sidebands - # sim.config.decompose_spectrum = "spin_system" - - # sim.spin_systems = spin_systems - # sim.methods = [method] - # sim.run(pack_as_csdm=False) - - # amp = sim.methods[0].simulation.real - print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - \ No newline at end of file diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index 32d92d3..c6bb444 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,7 +1,7 @@ -from mrinversion.kernel.csa_aniso import DAS # NOQA -from mrinversion.kernel.csa_aniso import MQMAS # NOQA -from mrinversion.kernel.csa_aniso import SL_MQMAS # NOQA -from mrinversion.kernel.csa_aniso import SL_MQMASnodist +from mrinversion.kernel.quad_aniso import DAS # NOQA +from mrinversion.kernel.quad_aniso import MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMASnodist from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 2f87bc0..3ec604b 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -1,3 +1,269 @@ +import warnings + +import csdmpy as cp +import numpy as np +from mrsimulator.utils import get_spectral_dimensions + +from mrsimulator import Simulator +from mrsimulator import SpinSystem +from mrsimulator.method import Method + +from mrinversion.kernel.base import LineShape + + +class DAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # "DAS", + ) + + def kernel(self, supersampling, eta_bound = 1): + # update method for DAS spectra events + das_event = dict( + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], + freq_contrib=["Quad2_4"], + ) + self.method_args["spectral_dimensions"][0]["events"] = [das_event] + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + return self._averaged_kernel(amp, supersampling) + + +class MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + + def kernel(self, supersampling, eta_bound = 1): + # update method for DAS spectra events + mqmas_events = [ + dict( + fraction= -9/50, + transition_queries=[{"ch1": {"P": [-3], "D": [0]}}] + ), + dict( + fraction= 27/50, + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}] + ), + ] + + self.method_args["spectral_dimensions"][0]["events"] = mqmas_events + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound, calc_pos=True) + + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + warnings.warn("This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative.") + + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMASnodist(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound = 1): + import sys + sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') + # import src.processing as smproc + import src.simulation as smsim + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray([smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type='c0_c4', + contribs='c0_c4', + return_array=True, + distorted=False) for mysys in spin_systems]) + + warnings.warn("This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative.") + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound = 1, cq_posneg=True): + import sys + sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') + # import src.processing as smproc + import src.simulation as smsim + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + + if eta_bound == 1 and cq_posneg: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound, cq_posneg) + + if eta_bound == 1: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) + for cq_, e,abun in zip(Cq, eta,abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray([smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type='c0_c4', + contribs='c0_c4', + return_array=True, + distorted=True) for mysys in spin_systems]) + + warnings.warn("This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative.") + return self._averaged_kernel(amp, supersampling) + + + + + # from copy import deepcopy # import numpy as np # from mrsimulator import Simulator diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index e19fe32..c610ac7 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -49,103 +49,44 @@ def x_y_to_zeta_eta(x, y): return zeta * x_unit, eta -def _x_y_to_zeta_eta(x, y, eta_bound=1): +def _x_y_to_zeta_eta(x, y, eta_bound=1,calc_pos=False): """Same as def x_y_to_zeta_eta, but for ndarrays.""" x = np.abs(x).ravel() y = np.abs(y).ravel() - # print(x/1e6) - # print(y/1e6) - n_lost = 0 zeta = np.sqrt(x**2 + y**2) # + offset abundances = np.ones(zeta.shape) eta = np.ones(zeta.shape) eta = eta.tolist() - # print(zeta/1e6) zeta = zeta.tolist() - # x.tolist() - # y.tolist() - - # print(f'eta length pre: {len(eta)}') - # print(eta.shape) del_these = [] for index,_ in enumerate(x): if x[index] > y[index]: - zeta[index] = - zeta[index] - this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) - if this_eta < eta_bound: - eta[index] = this_eta + if not calc_pos: + zeta[index] = - zeta[index] + this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) + if this_eta < eta_bound: + eta[index] = this_eta + else: + abundances[index] = 0 else: - # del_these.append(index) - # del zeta[index], eta[index] - # n_lost += 1 abundances[index] = 0 + elif x[index] < y[index]: this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) if this_eta < eta_bound: eta[index] = this_eta else: - # del zeta[index], eta[index] - # del_these.append(index) - # n_lost += 1 abundances[index] = 0 elif x[index] == y[index] and eta_bound < 1: - # del_these.append(index) - # n_lost += 1 abundances[index] = 0 - # print(np.asarray(eta)) - # print(del_these) - # print(f'eta length post: {len(eta)}') - # print(f'n lost: {n_lost}') - if eta_bound ==1: + if eta_bound ==1 and not calc_pos: return np.asarray(zeta), np.asarray(eta) else: for idx in del_these[::-1]: del zeta[idx], eta[idx] - # print(np.asarray(eta)) return np.asarray(zeta), np.asarray(eta), abundances - - - - -# def _x_y_to_zeta_eta(x, y, eta_bound=1): -# """Same as def x_y_to_zeta_eta, but for ndarrays.""" -# x = np.abs(x) -# y = np.abs(y) -# # print(x) -# # print(y) -# n_lost = 0 -# zeta = np.sqrt(x**2 + y**2) # + offset -# # print(zeta) -# eta = np.ones(zeta.shape)*100 -# # print(eta.shape) -# for index,_ in np.ndenumerate(x): -# if x[index] > y[index]: -# zeta[index] = - zeta[index] -# this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) -# if this_eta < eta_bound: -# eta[index] = this_eta -# else: -# n_lost += 1 -# elif x[index] < y[index]: -# this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) -# if this_eta < eta_bound: -# eta[index] = this_eta -# else: -# n_lost += 1 -# print(eta) - -# if eta_bound ==1: -# return zeta.ravel(), eta.ravel() -# else: -# return zeta.ravel(), eta.ravel(), n_lost - - - - - - def zeta_eta_to_x_y(zeta, eta): r"""Convert the coordinates :math:`(\zeta,\eta)` to :math:`(x, y)` using the following definition, @@ -190,7 +131,7 @@ def zeta_eta_to_x_y(zeta, eta): return x.ravel(), y.ravel() -def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound): +def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound,calc_pos=False): """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) @@ -207,27 +148,27 @@ def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound): np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" ) - return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound) + return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound, calc_pos) -def _x_y_to_cq_eta_distribution(grid, supersampling): - """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" - x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) - y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) +# def _x_y_to_cq_eta_distribution(grid, supersampling): +# """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" +# x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) +# y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) - if x_coordinates.unit.physical_type == "frequency": - x_coordinates = x_coordinates.to("Hz").value - y_coordinates = y_coordinates.to("Hz").value +# if x_coordinates.unit.physical_type == "frequency": +# x_coordinates = x_coordinates.to("Hz").value +# y_coordinates = y_coordinates.to("Hz").value - elif x_coordinates.unit.physical_type == "dimensionless": - x_coordinates = x_coordinates.to("ppm").value - y_coordinates = y_coordinates.to("ppm").value +# elif x_coordinates.unit.physical_type == "dimensionless": +# x_coordinates = x_coordinates.to("ppm").value +# y_coordinates = y_coordinates.to("ppm").value - x_mesh, y_mesh = np.meshgrid( - np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" - ) +# x_mesh, y_mesh = np.meshgrid( +# np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" +# ) - return _x_y_to_cq_eta(x_mesh, y_mesh) +# return _x_y_to_cq_eta(x_mesh, y_mesh) def _supersampled_coordinates(dimension, supersampling=1): @@ -274,52 +215,47 @@ def _supersampled_coordinates(dimension, supersampling=1): return array -def cq_eta_to_x_y(cq, eta): - r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the - following definition, +# def cq_eta_to_x_y(cq, eta): +# r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the +# following definition, - .. math:: - \begin{array}{rl} - x &= C_q \sin\theta, \\ - y &= C_q \cos\theta - \end{array} {~~~~~~~~} +# .. math:: +# \begin{array}{rl} +# x &= C_q \sin\theta, \\ +# y &= C_q \cos\theta +# \end{array} {~~~~~~~~} - where :math:`\theta = \frac{\pi}{2}\eta`. +# where :math:`\theta = \frac{\pi}{2}\eta`. - Args: - x: ndarray or list of floats. The coordinate x. - y: ndarray or list of floats. The coordinate y. +# Args: +# x: ndarray or list of floats. The coordinate x. +# y: ndarray or list of floats. The coordinate y. - Return: - A list of ndarrays. The first array holds the coordinate :math:`x`. The - second array holds the coordinates :math:`y`. - """ - cq = np.asarray(cq) - eta = np.asarray(eta) +# Return: +# A list of ndarrays. The first array holds the coordinate :math:`x`. The +# second array holds the coordinates :math:`y`. +# """ +# cq = np.asarray(cq) +# eta = np.asarray(eta) - theta = np.pi * eta / 2.0 - x = np.zeros(cq.size) - y = np.zeros(cq.size) +# theta = np.pi * eta / 2.0 +# x = np.zeros(cq.size) +# y = np.zeros(cq.size) - index = np.arange(len(cq)) - x[index] = cq[index] * np.cos(theta[index]) - y[index] = cq[index] * np.sin(theta[index]) +# index = np.arange(len(cq)) +# x[index] = cq[index] * np.cos(theta[index]) +# y[index] = cq[index] * np.sin(theta[index]) - return x.ravel(), y.ravel() +# return x.ravel(), y.ravel() -def _x_y_to_cq_eta(x, y): - """Same as def x_y_to_zeta_eta, but for ndarrays.""" - x = np.abs(x) - y = np.abs(y) - cq = np.sqrt(x**2 + y**2) # + offset - eta = np.ones(cq.shape) - # index = np.where(x > y) - index = np.arange(len(cq)) - # zeta[index] = -zeta[index] - eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) - - # index = np.where(x < y) - # eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index]) - - return cq.ravel(), eta.ravel() \ No newline at end of file +# def _x_y_to_cq_eta(x, y): +# """Same as def x_y_to_zeta_eta, but for ndarrays.""" +# x = np.abs(x) +# y = np.abs(y) +# cq = np.sqrt(x**2 + y**2) # + offset +# eta = np.ones(cq.shape) +# index = np.arange(len(cq)) +# eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) + +# return cq.ravel(), eta.ravel() \ No newline at end of file diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 81b7458..6ca1aa7 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -60,7 +60,11 @@ def __init__( regularizer=None, inverse_dimension=None, method="gradient_decent", - ): + xygrid='full' + ): + + if xygrid != 'full' and xygrid != 'mirrored' and xygrid != 'positive' and xygrid != 'negative': + raise Exception("Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'.") self.hyperparameters = {"lambda": lambda1, "alpha": alpha} self.max_iterations = max_iterations @@ -74,6 +78,7 @@ def __init__( # attributes self.f = None self.n_iter = None + self.xygrid = xygrid def fit(self, K, s): r"""Fit the model using the coordinate descent method from scikit-learn. @@ -194,6 +199,45 @@ def _sol_to_csdm(self, s): self.f.dimensions[1] = s.dimensions[1] self.f.dimensions[0] = self.inverse_dimension[0] + if self.xygrid != 'full': + if self.xygrid == 'mirrored': + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == 'negative': + self.f = self.flip_inversion() + + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx,val in np.ndenumerate(y): + dim0=idx[0] + dim1=idx[1] + if dim1 < dim0: + if y[dim1,dim0] < 1e-6: + y[dim1,dim0] = val + y[dim0,dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, + dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx,val in np.ndenumerate(y): + dim0=idx[0] + dim1=idx[1] + if dim1 < dim0: + if y[dim1,dim0] < 1e-6: + y[dim1,dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, + dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + def predict(self, K): r"""Predict the signal using the linear model. @@ -271,7 +315,10 @@ def __init__( inverse_dimension=None, n_jobs=-1, method="gradient_decent", - ): + xygrid='full' + ): + if xygrid != 'full' and xygrid != 'mirrored' and xygrid != 'positive' and xygrid != 'negative': + raise Exception("Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'.") self.cv_alphas = np.asarray(alphas).ravel() self.cv_lambdas = np.asarray(lambdas).ravel() @@ -292,8 +339,10 @@ def __init__( self.verbose = verbose self.inverse_dimension = inverse_dimension self.f_shape = tuple(item.count for item in inverse_dimension)[::-1] + self.xygrid = xygrid - def fit(self, K, s): + + def fit(self, K, s, cv_map_as_csdm=True): r"""Fit the model using the coordinate descent method from scikit-learn for all alpha anf lambda values using the `n`-folds cross-validation technique. The cross-validation metric is the mean squared error. @@ -383,19 +432,46 @@ def fit(self, K, s): # Calculate the solution using the complete data at the optimized lambda and # alpha values - self.opt = GeneralL2Lasso( - alpha=self.hyperparameters["alpha"], - lambda1=self.hyperparameters["lambda"], - max_iterations=self.max_iterations, - tolerance=self.tolerance, - positive=self.positive, - regularizer=self.regularizer, - inverse_dimension=self.inverse_dimension, - method=self.method, - ) - self.opt.fit(K, s) - self.f = self.opt.f - self.cv_map_to_csdm() + if self.xygrid == 'full': + self.opt = GeneralL2Lasso( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + ) + self.opt.fit(K, s) + self.f = self.opt.f + if cv_map_as_csdm: + self.cv_map_to_csdm() + + else: + self.opt = GeneralL2Lasso( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + xygrid=self.xygrid + ) + + self.opt.fit(K, s) + self.f = self.opt.f + + if self.xygrid == 'mirrored': + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == 'negative': + self.f = self.flip_inversion() + + if cv_map_as_csdm: + self.cv_map_to_csdm() def cv_map_to_csdm(self): # convert cv_map to csdm @@ -476,6 +552,37 @@ def _get_minimizer(self): random_state=None, ) + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx,val in np.ndenumerate(y): + dim0=idx[0] + dim1=idx[1] + if dim1 < dim0: + if y[dim1,dim0] < 1e-6: + y[dim1,dim0] = val + y[dim0,dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, + dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx,val in np.ndenumerate(y): + dim0=idx[0] + dim1=idx[1] + if dim1 < dim0: + if y[dim1,dim0] < 1e-6: + y[dim1,dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, + dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + def predict(self, K): r"""Predict the signal using the linear model. diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index 9f9b202..3c7f5ed 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -1,6 +1,8 @@ from ._base_l1l2 import GeneralL2Lasso from ._base_l1l2 import GeneralL2LassoCV +import warnings + __author__ = "Deepansh J. Srivastava" __email__ = "srivastava.89@osu.edu" @@ -59,6 +61,15 @@ class SmoothLasso(GeneralL2Lasso): If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only positive values, else the solution may contain positive and negative amplitudes. The default is True. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. Attributes ---------- @@ -80,6 +91,7 @@ def __init__( tolerance=1e-5, positive=True, method="gradient_decent", + xygrid='full' ): super().__init__( alpha=alpha, @@ -90,6 +102,7 @@ def __init__( regularizer="smooth lasso", inverse_dimension=inverse_dimension, method=method, + xygrid=xygrid ) @@ -168,6 +181,15 @@ class SmoothLassoCV(GeneralL2LassoCV): n_jobs: int The number of CPUs used for computation. The default is -1, that is, all available CPUs are used. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. Attributes @@ -199,6 +221,7 @@ def __init__( verbose=False, n_jobs=-1, method="gradient_decent", + xygrid='full' ): super().__init__( alphas=alphas, @@ -215,4 +238,5 @@ def __init__( verbose=verbose, n_jobs=n_jobs, method=method, + xygrid=xygrid ) diff --git a/mrinversion/linear_model/utils.py b/mrinversion/linear_model/utils.py new file mode 100644 index 0000000..f0d17ea --- /dev/null +++ b/mrinversion/linear_model/utils.py @@ -0,0 +1,199 @@ +from scipy.optimize import minimize +import numpy as np +from mrinversion.linear_model import SmoothLassoCV + +class CVMinimizer: + def __init__( + self, + inverse_dimension, + compressed_K, + compressed_s, + guess_neglogalpha, + guess_negloglambda, + sigma=0.0, + folds=10, + xygrid='full' + ): + self.inverse_dimension = inverse_dimension + self.compressed_K = compressed_K + self.compressed_s = compressed_s + self.guess = np.array([guess_neglogalpha, guess_negloglambda]) + self.sigma = sigma + self.folds = folds + self.xygrid = xygrid + + @staticmethod + def _cv_min_withprint(loghyperparameters,self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + print(f'alpha: {-np.log10(alpha)}, lambda: {-np.log10(lam)}, cv: {cv}') + return cv + + @staticmethod + def _cv_min_noprint(loghyperparameters, self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + return cv + + @staticmethod + def _print_result(this_result, print_steps, prev_result=None): + if prev_result: + if this_result.fun < prev_result.fun: + if this_result.success: + print('This minimization improved on the previous result and completed the minimization in the iterations given. ') + else: + print('This minimization improved on the previous result, but did not finish in the given number of iterations.') + else: + if this_result.success: + print('This minimization minimized to a worse solution.') + else: + print('This minimization did not improve on the previous result, and did not finish in the given number of iterations.') + if print_steps: + print(f'This minimization results:\n\n {this_result}') + else: + if this_result.success: + print('This minimization completed in the iterations given. ') + else: + print('This minimization did not finish in the given number of iterations.') + if print_steps: + print(f'This minimization results:\n\n {this_result}') + + print('--------------------------------------------------\n') + + def _minimize(self, guess, tol, maxiter, print_funcevals, print_steps, prev_result=None): + if print_funcevals: + simplex_res = minimize( + fun=self._cv_min_withprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={'maxiter': maxiter}, + method='Nelder-Mead' + ) + else: + simplex_res = minimize( + fun=self._cv_min_noprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={'maxiter': maxiter}, + method='Nelder-Mead' + ) + self._print_result(this_result=simplex_res, print_steps=print_steps, prev_result=prev_result) + return simplex_res + + def mincv_close(self, tol=0.0002, maxiter=500, print_steps=True, print_funcevals=False): + simplex_res1 = self._minimize( + guess=self.guess, + tol=tol*1.5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps + ) + + if simplex_res1.success: + simplex_res2 = self._minimize( + guess=simplex_res1.x, + tol=tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + prev_result=simplex_res1 + ) + if simplex_res2.fun < simplex_res1.fun: + return simplex_res2 + else: + return simplex_res1 + else: + return simplex_res1 + + def mincv(self, tol=0.0002, maxiter=100, print_steps=True, print_funcevals=False, guess_bounds=[[2,10],[2,10]], guesses='all'): + ''' options for guesses: "all", "cartesian", "diagonals"''' + lowval_alpha = np.average([self.guess[0], guess_bounds[0][0]]) + highval_alpha = np.average([self.guess[0], guess_bounds[0][1]]) + lowval_lambda = np.average([self.guess[1], guess_bounds[1][0]]) + highval_lambda = np.average([self.guess[1], guess_bounds[1][1]]) + if guesses == 'all': + trythese = np.zeros((9,2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, highval_lambda] + trythese[2] = [highval_alpha, self.guess[1]] + trythese[3] = [highval_alpha, lowval_lambda] + trythese[4] = [self.guess[0], lowval_lambda] + trythese[5] = [lowval_alpha, lowval_lambda] + trythese[6] = [lowval_alpha, self.guess[1]] + trythese[7] = [lowval_alpha, highval_lambda] + trythese[8] = self.guess + elif guesses == 'cartesian': + trythese = np.zeros((5,2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, self.guess[1]] + trythese[2] = [self.guess[0], lowval_lambda] + trythese[3] = [lowval_alpha, self.guess[1]] + trythese[4] = self.guess + elif guesses == 'diagonal': + trythese = np.zeros((5,2)) + trythese[0] = [highval_alpha, highval_lambda] + trythese[1] = [highval_alpha, lowval_lambda] + trythese[2] = [lowval_alpha, lowval_lambda] + trythese[3] = [lowval_alpha, highval_lambda] + trythese[4] = self.guess + else: + print('choose a valid choice for "guesses"') + return None + print(f'Starting points to try: {trythese}') + first_step = [ + self._minimize( + guess=thisguess, + tol = tol*5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps + ) + for thisguess in trythese + ] + print('-----------------------------------------------------------------') + mins = np.asarray([step.fun for step in first_step]) + min_idx = np.where(mins == mins.min())[0][0] + + final_step = self._minimize( + guess = first_step[min_idx].x, + tol = tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=False + ) + + self._print_result(final_step, print_steps=True, prev_result=first_step[min_idx]) + + if final_step.fun > first_step[min_idx].fun: + print('Second minimization did not improve.') + return first_step[min_idx] + else: + return final_step From 94c3f94ac750c1312bdcd45858d2d88869f9b611 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:28:39 -0400 Subject: [PATCH 08/20] add CV minimization and xygrid parameter to SmoothLasso and SmoothLassoCV --- mrinversion/kernel/base.py | 13 +- mrinversion/kernel/csa_aniso.py | 297 +--------------------- mrinversion/kernel/nmr.py | 10 +- mrinversion/kernel/quad_aniso.py | 300 +++++++++++++++++++++++ mrinversion/kernel/utils.py | 187 +++++--------- mrinversion/linear_model/_base_l1l2.py | 145 +++++++++-- mrinversion/linear_model/smooth_lasso.py | 24 ++ mrinversion/linear_model/utils.py | 225 +++++++++++++++++ 8 files changed, 755 insertions(+), 446 deletions(-) create mode 100644 mrinversion/linear_model/utils.py diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 9c48b2b..d21c71a 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -89,7 +89,7 @@ def _averaged_kernel(self, amp, supersampling, xy_grid=True, mask_kernel=False): # for i in range(len(K)): # mask[i][len(K[0])-i:] += 1 # K = ma.masked_array(K, mask=mask) - + inv_size = np.asarray([item.count for item in inverse_kernel_dimension]).prod() K = K.reshape(inv_size, self.kernel_dimension.count).T return K @@ -165,17 +165,16 @@ def __init__( if number_of_sidebands is None: self.number_of_sidebands = dim.count - def _get_zeta_eta(self, supersampling, eta_bound=1): + def _get_zeta_eta(self, supersampling, eta_bound=1, calc_pos=False): """Return zeta and eta coordinates over x-y grid""" - - if eta_bound == 1: + if eta_bound == 1 and not calc_pos: zeta, eta = _x_y_to_zeta_eta_distribution( self.inverse_kernel_dimension, supersampling, eta_bound ) return zeta, eta else: - zeta, eta, abundances= _x_y_to_zeta_eta_distribution( - self.inverse_kernel_dimension, supersampling, eta_bound + zeta, eta, abundances = _x_y_to_zeta_eta_distribution( + self.inverse_kernel_dimension, supersampling, eta_bound, calc_pos ) return zeta, eta, abundances @@ -226,4 +225,4 @@ def _check_dimension_type(dimensions, direction, dimension_quantity, kernel_type f"is required for the `{kernel_type}` kernel, instead got " f"`{item.quantity_name}` as the quantity name for the dimension at " f"index {i}." - ) \ No newline at end of file + ) diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index 5a391da..29ac667 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -2,12 +2,12 @@ import csdmpy as cp import numpy as np -from mrsimulator.utils import get_spectral_dimensions - from mrsimulator import Simulator from mrsimulator import SpinSystem -from mrsimulator.method import Method, SpectralEvent +from mrsimulator.method import Method +from mrsimulator.method import SpectralEvent from mrsimulator.method.lib import BlochDecaySpectrum +from mrsimulator.utils import get_spectral_dimensions from mrinversion.kernel.base import LineShape @@ -184,294 +184,3 @@ def __init__( None, None, ) - - -class DAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # "DAS", - ) - - def kernel(self, supersampling, eta_bound = 1): - # update method for DAS spectra events - das_event = dict( - transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], - freq_contrib=["Quad2_4"], - ) - self.method_args["spectral_dimensions"][0]["events"] = [das_event] - - method = Method.parse_dict_with_units(self.method_args) - isotope = self.method_args["channels"][0] - - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - sim = Simulator() - sim.config.number_of_sidebands = self.number_of_sidebands - sim.config.decompose_spectrum = "spin_system" - - sim.spin_systems = spin_systems - sim.methods = [method] - sim.run(pack_as_csdm=False) - - amp = sim.methods[0].simulation.real - - return self._averaged_kernel(amp, supersampling) - - -class MQMAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # "DAS", - ) - - def kernel(self, supersampling, eta_bound = 1): - # update method for DAS spectra events - mqmas_events = [ - dict( - fraction= -9/50, - transition_queries=[{"ch1": {"P": [-3], "D": [0]}}] - ), - dict( - fraction= 27/50, - transition_queries=[{"ch1": {"P": [-1], "D": [0]}}] - ), - ] - - self.method_args["spectral_dimensions"][0]["events"] = mqmas_events - - method = Method.parse_dict_with_units(self.method_args) - isotope = self.method_args["channels"][0] - - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - sim = Simulator() - sim.config.number_of_sidebands = self.number_of_sidebands - sim.config.decompose_spectrum = "spin_system" - - sim.spin_systems = spin_systems - sim.methods = [method] - sim.run(pack_as_csdm=False) - - amp = sim.methods[0].simulation.real - print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - - -class SL_MQMASnodist(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - exp_dict, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # - # "DAS", - ) - self.exp_dict = exp_dict - self.anisotropic_dimension = anisotropic_dimension - - def kernel(self, supersampling, eta_bound = 1): - import sys - sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') - # import src.processing as smproc - import src.simulation as smsim - # import src.fitting as smfit - - isotope = self.method_args["channels"][0] - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - # print(f'Cq: {Cq}') - # print(f'eta: {eta}') - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - # print(self.anisotropic_dimension) - - obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) - spec_dim = get_spectral_dimensions(obj) - # print(obj.x[0]) - # print(self.anisotropic_dimension) - # print(spec_dim) - - amp = np.asarray([smsim.simulate_onesite_lineshape( - self.exp_dict, - mysys, - spec_dim[0], - input_type='c0_c4', - contribs='c0_c4', - return_array=True, - distorted=False) for mysys in spin_systems]) - # sim = Simulator() - # sim.config.number_of_sidebands = self.number_of_sidebands - # sim.config.decompose_spectrum = "spin_system" - - # sim.spin_systems = spin_systems - # sim.methods = [method] - # sim.run(pack_as_csdm=False) - # obj_pre = cp.CSDM( - # dimensions=[ - # aniso - # ] - # ) - # amp = sim.methods[0].simulation.real - # print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - - -class SL_MQMAS(LineShape): - def __init__( - self, - anisotropic_dimension, - inverse_kernel_dimension, - channel, - exp_dict, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="600 Hz", - number_of_sidebands=None, - ): - super().__init__( - anisotropic_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - # - # "DAS", - ) - self.exp_dict = exp_dict - self.anisotropic_dimension = anisotropic_dimension - - def kernel(self, supersampling, eta_bound = 1): - import sys - sys.path.insert(0, '/home/lexicon2810/github-repos-WSL/mrsmqmas') - # import src.processing as smproc - import src.simulation as smsim - # import src.fitting as smfit - - isotope = self.method_args["channels"][0] - if eta_bound == 1: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) - - print(f'Cq: {Cq}') - print(f'eta: {eta}') - if eta_bound == 1: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], abundance=abun) - for cq_, e,abun in zip(Cq, eta,abundances) - ] - # print(self.anisotropic_dimension) - - obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) - spec_dim = get_spectral_dimensions(obj) - - amp = np.asarray([smsim.simulate_onesite_lineshape( - self.exp_dict, - mysys, - spec_dim[0], - input_type='c0_c4', - contribs='c0_c4', - return_array=True, - distorted=True) for mysys in spin_systems]) - # sim = Simulator() - # sim.config.number_of_sidebands = self.number_of_sidebands - # sim.config.decompose_spectrum = "spin_system" - - # sim.spin_systems = spin_systems - # sim.methods = [method] - # sim.run(pack_as_csdm=False) - - # amp = sim.methods[0].simulation.real - print(f'amp shape: {amp.shape}') - return self._averaged_kernel(amp, supersampling) - \ No newline at end of file diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index 32d92d3..b2206aa 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,9 +1,9 @@ -from mrinversion.kernel.csa_aniso import DAS # NOQA -from mrinversion.kernel.csa_aniso import MQMAS # NOQA -from mrinversion.kernel.csa_aniso import SL_MQMAS # NOQA -from mrinversion.kernel.csa_aniso import SL_MQMASnodist from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA +from mrinversion.kernel.quad_aniso import DAS # NOQA +from mrinversion.kernel.quad_aniso import MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMASnodist -# from mrinversion.kernel.quad_aniso import MQMAS # NOQA \ No newline at end of file +# from mrinversion.kernel.quad_aniso import MQMAS # NOQA diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 2f87bc0..9a40e3d 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -1,3 +1,303 @@ +import warnings + +import csdmpy as cp +import numpy as np +from mrsimulator import Simulator +from mrsimulator import SpinSystem +from mrsimulator.method import Method +from mrsimulator.utils import get_spectral_dimensions + +from mrinversion.kernel.base import LineShape + + +class DAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # "DAS", + ) + + def kernel(self, supersampling, eta_bound=1): + # update method for DAS spectra events + das_event = dict( + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], + freq_contrib=["Quad2_4"], + ) + self.method_args["spectral_dimensions"][0]["events"] = [das_event] + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] + ) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + return self._averaged_kernel(amp, supersampling) + + +class MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + + def kernel(self, supersampling, eta_bound=1): + # update method for DAS spectra events + mqmas_events = [ + dict(fraction=-9 / 50, transition_queries=[{"ch1": {"P": [-3], "D": [0]}}]), + dict(fraction=27 / 50, transition_queries=[{"ch1": {"P": [-1], "D": [0]}}]), + ] + + self.method_args["spectral_dimensions"][0]["events"] = mqmas_events + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + Cq, eta, abundances = self._get_zeta_eta( + supersampling, eta_bound, calc_pos=True + ) + + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + ) + + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMASnodist(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound=1): + import sys + + sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") + # import src.processing as smproc + import src.simulation as smsim + + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] + ) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray( + [ + smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type="c0_c4", + contribs="c0_c4", + return_array=True, + distorted=False, + ) + for mysys in spin_systems + ] + ) + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + ) + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound=1, cq_posneg=True): + import sys + + sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") + # import src.processing as smproc + import src.simulation as smsim + + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + + if eta_bound == 1 and cq_posneg: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta( + supersampling, eta_bound, cq_posneg + ) + + if eta_bound == 1: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] + ) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray( + [ + smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type="c0_c4", + contribs="c0_c4", + return_array=True, + distorted=True, + ) + for mysys in spin_systems + ] + ) + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + ) + return self._averaged_kernel(amp, supersampling) + + # from copy import deepcopy # import numpy as np # from mrsimulator import Simulator diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index e19fe32..83ae64a 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -49,103 +49,45 @@ def x_y_to_zeta_eta(x, y): return zeta * x_unit, eta -def _x_y_to_zeta_eta(x, y, eta_bound=1): +def _x_y_to_zeta_eta(x, y, eta_bound=1, calc_pos=False): """Same as def x_y_to_zeta_eta, but for ndarrays.""" x = np.abs(x).ravel() y = np.abs(y).ravel() - # print(x/1e6) - # print(y/1e6) - n_lost = 0 zeta = np.sqrt(x**2 + y**2) # + offset abundances = np.ones(zeta.shape) eta = np.ones(zeta.shape) eta = eta.tolist() - # print(zeta/1e6) zeta = zeta.tolist() - # x.tolist() - # y.tolist() - - # print(f'eta length pre: {len(eta)}') - # print(eta.shape) del_these = [] - for index,_ in enumerate(x): + for index, _ in enumerate(x): if x[index] > y[index]: - zeta[index] = - zeta[index] - this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) - if this_eta < eta_bound: + if not calc_pos: + zeta[index] = -zeta[index] + this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) + if this_eta < eta_bound: eta[index] = this_eta + else: + abundances[index] = 0 else: - # del_these.append(index) - # del zeta[index], eta[index] - # n_lost += 1 abundances[index] = 0 + elif x[index] < y[index]: this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) if this_eta < eta_bound: eta[index] = this_eta else: - # del zeta[index], eta[index] - # del_these.append(index) - # n_lost += 1 abundances[index] = 0 elif x[index] == y[index] and eta_bound < 1: - # del_these.append(index) - # n_lost += 1 abundances[index] = 0 - - # print(np.asarray(eta)) - # print(del_these) - # print(f'eta length post: {len(eta)}') - # print(f'n lost: {n_lost}') - if eta_bound ==1: + + if eta_bound == 1 and not calc_pos: return np.asarray(zeta), np.asarray(eta) else: for idx in del_these[::-1]: del zeta[idx], eta[idx] - # print(np.asarray(eta)) return np.asarray(zeta), np.asarray(eta), abundances - - - -# def _x_y_to_zeta_eta(x, y, eta_bound=1): -# """Same as def x_y_to_zeta_eta, but for ndarrays.""" -# x = np.abs(x) -# y = np.abs(y) -# # print(x) -# # print(y) -# n_lost = 0 -# zeta = np.sqrt(x**2 + y**2) # + offset -# # print(zeta) -# eta = np.ones(zeta.shape)*100 -# # print(eta.shape) -# for index,_ in np.ndenumerate(x): -# if x[index] > y[index]: -# zeta[index] = - zeta[index] -# this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) -# if this_eta < eta_bound: -# eta[index] = this_eta -# else: -# n_lost += 1 -# elif x[index] < y[index]: -# this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) -# if this_eta < eta_bound: -# eta[index] = this_eta -# else: -# n_lost += 1 -# print(eta) - -# if eta_bound ==1: -# return zeta.ravel(), eta.ravel() -# else: -# return zeta.ravel(), eta.ravel(), n_lost - - - - - - def zeta_eta_to_x_y(zeta, eta): r"""Convert the coordinates :math:`(\zeta,\eta)` to :math:`(x, y)` using the following definition, @@ -190,7 +132,7 @@ def zeta_eta_to_x_y(zeta, eta): return x.ravel(), y.ravel() -def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound): +def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound, calc_pos=False): """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) @@ -207,27 +149,27 @@ def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound): np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" ) - return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound) + return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound, calc_pos) -def _x_y_to_cq_eta_distribution(grid, supersampling): - """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" - x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) - y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) +# def _x_y_to_cq_eta_distribution(grid, supersampling): +# """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" +# x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) +# y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) - if x_coordinates.unit.physical_type == "frequency": - x_coordinates = x_coordinates.to("Hz").value - y_coordinates = y_coordinates.to("Hz").value +# if x_coordinates.unit.physical_type == "frequency": +# x_coordinates = x_coordinates.to("Hz").value +# y_coordinates = y_coordinates.to("Hz").value - elif x_coordinates.unit.physical_type == "dimensionless": - x_coordinates = x_coordinates.to("ppm").value - y_coordinates = y_coordinates.to("ppm").value +# elif x_coordinates.unit.physical_type == "dimensionless": +# x_coordinates = x_coordinates.to("ppm").value +# y_coordinates = y_coordinates.to("ppm").value - x_mesh, y_mesh = np.meshgrid( - np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" - ) +# x_mesh, y_mesh = np.meshgrid( +# np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" +# ) - return _x_y_to_cq_eta(x_mesh, y_mesh) +# return _x_y_to_cq_eta(x_mesh, y_mesh) def _supersampled_coordinates(dimension, supersampling=1): @@ -274,52 +216,47 @@ def _supersampled_coordinates(dimension, supersampling=1): return array -def cq_eta_to_x_y(cq, eta): - r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the - following definition, +# def cq_eta_to_x_y(cq, eta): +# r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the +# following definition, - .. math:: - \begin{array}{rl} - x &= C_q \sin\theta, \\ - y &= C_q \cos\theta - \end{array} {~~~~~~~~} +# .. math:: +# \begin{array}{rl} +# x &= C_q \sin\theta, \\ +# y &= C_q \cos\theta +# \end{array} {~~~~~~~~} - where :math:`\theta = \frac{\pi}{2}\eta`. +# where :math:`\theta = \frac{\pi}{2}\eta`. - Args: - x: ndarray or list of floats. The coordinate x. - y: ndarray or list of floats. The coordinate y. +# Args: +# x: ndarray or list of floats. The coordinate x. +# y: ndarray or list of floats. The coordinate y. - Return: - A list of ndarrays. The first array holds the coordinate :math:`x`. The - second array holds the coordinates :math:`y`. - """ - cq = np.asarray(cq) - eta = np.asarray(eta) +# Return: +# A list of ndarrays. The first array holds the coordinate :math:`x`. The +# second array holds the coordinates :math:`y`. +# """ +# cq = np.asarray(cq) +# eta = np.asarray(eta) - theta = np.pi * eta / 2.0 - x = np.zeros(cq.size) - y = np.zeros(cq.size) +# theta = np.pi * eta / 2.0 +# x = np.zeros(cq.size) +# y = np.zeros(cq.size) - index = np.arange(len(cq)) - x[index] = cq[index] * np.cos(theta[index]) - y[index] = cq[index] * np.sin(theta[index]) +# index = np.arange(len(cq)) +# x[index] = cq[index] * np.cos(theta[index]) +# y[index] = cq[index] * np.sin(theta[index]) - return x.ravel(), y.ravel() +# return x.ravel(), y.ravel() -def _x_y_to_cq_eta(x, y): - """Same as def x_y_to_zeta_eta, but for ndarrays.""" - x = np.abs(x) - y = np.abs(y) - cq = np.sqrt(x**2 + y**2) # + offset - eta = np.ones(cq.shape) - # index = np.where(x > y) - index = np.arange(len(cq)) - # zeta[index] = -zeta[index] - eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) - - # index = np.where(x < y) - # eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index]) - - return cq.ravel(), eta.ravel() \ No newline at end of file +# def _x_y_to_cq_eta(x, y): +# """Same as def x_y_to_zeta_eta, but for ndarrays.""" +# x = np.abs(x) +# y = np.abs(y) +# cq = np.sqrt(x**2 + y**2) # + offset +# eta = np.ones(cq.shape) +# index = np.arange(len(cq)) +# eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) + +# return cq.ravel(), eta.ravel() diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 81b7458..256012c 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -60,8 +60,19 @@ def __init__( regularizer=None, inverse_dimension=None, method="gradient_decent", + xygrid="full", ): + if ( + xygrid != "full" + and xygrid != "mirrored" + and xygrid != "positive" + and xygrid != "negative" + ): + raise Exception( + "Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'." + ) + self.hyperparameters = {"lambda": lambda1, "alpha": alpha} self.max_iterations = max_iterations self.tolerance = tolerance @@ -74,6 +85,7 @@ def __init__( # attributes self.f = None self.n_iter = None + self.xygrid = xygrid def fit(self, K, s): r"""Fit the model using the coordinate descent method from scikit-learn. @@ -194,6 +206,42 @@ def _sol_to_csdm(self, s): self.f.dimensions[1] = s.dimensions[1] self.f.dimensions[0] = self.inverse_dimension[0] + if self.xygrid != "full": + if self.xygrid == "mirrored": + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() + + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + y[dim0, dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + def predict(self, K): r"""Predict the signal using the linear model. @@ -271,7 +319,17 @@ def __init__( inverse_dimension=None, n_jobs=-1, method="gradient_decent", + xygrid="full", ): + if ( + xygrid != "full" + and xygrid != "mirrored" + and xygrid != "positive" + and xygrid != "negative" + ): + raise Exception( + "Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'." + ) self.cv_alphas = np.asarray(alphas).ravel() self.cv_lambdas = np.asarray(lambdas).ravel() @@ -292,8 +350,9 @@ def __init__( self.verbose = verbose self.inverse_dimension = inverse_dimension self.f_shape = tuple(item.count for item in inverse_dimension)[::-1] + self.xygrid = xygrid - def fit(self, K, s): + def fit(self, K, s, cv_map_as_csdm=True): r"""Fit the model using the coordinate descent method from scikit-learn for all alpha anf lambda values using the `n`-folds cross-validation technique. The cross-validation metric is the mean squared error. @@ -383,19 +442,46 @@ def fit(self, K, s): # Calculate the solution using the complete data at the optimized lambda and # alpha values - self.opt = GeneralL2Lasso( - alpha=self.hyperparameters["alpha"], - lambda1=self.hyperparameters["lambda"], - max_iterations=self.max_iterations, - tolerance=self.tolerance, - positive=self.positive, - regularizer=self.regularizer, - inverse_dimension=self.inverse_dimension, - method=self.method, - ) - self.opt.fit(K, s) - self.f = self.opt.f - self.cv_map_to_csdm() + if self.xygrid == "full": + self.opt = GeneralL2Lasso( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + ) + self.opt.fit(K, s) + self.f = self.opt.f + if cv_map_as_csdm: + self.cv_map_to_csdm() + + else: + self.opt = GeneralL2Lasso( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + xygrid=self.xygrid, + ) + + self.opt.fit(K, s) + self.f = self.opt.f + + if self.xygrid == "mirrored": + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() + + if cv_map_as_csdm: + self.cv_map_to_csdm() def cv_map_to_csdm(self): # convert cv_map to csdm @@ -476,6 +562,35 @@ def _get_minimizer(self): random_state=None, ) + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + y[dim0, dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + def predict(self, K): r"""Predict the signal using the linear model. @@ -674,4 +789,4 @@ def prepare_signal(sig): scale = np.sqrt(np.mean(np.abs(s_) ** 2)) s_ = s_ / scale - return s_, scale \ No newline at end of file + return s_, scale diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index 9f9b202..fcc13ef 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -1,3 +1,5 @@ +import warnings + from ._base_l1l2 import GeneralL2Lasso from ._base_l1l2 import GeneralL2LassoCV @@ -59,6 +61,15 @@ class SmoothLasso(GeneralL2Lasso): If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only positive values, else the solution may contain positive and negative amplitudes. The default is True. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. Attributes ---------- @@ -80,6 +91,7 @@ def __init__( tolerance=1e-5, positive=True, method="gradient_decent", + xygrid="full", ): super().__init__( alpha=alpha, @@ -90,6 +102,7 @@ def __init__( regularizer="smooth lasso", inverse_dimension=inverse_dimension, method=method, + xygrid=xygrid, ) @@ -168,6 +181,15 @@ class SmoothLassoCV(GeneralL2LassoCV): n_jobs: int The number of CPUs used for computation. The default is -1, that is, all available CPUs are used. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. Attributes @@ -199,6 +221,7 @@ def __init__( verbose=False, n_jobs=-1, method="gradient_decent", + xygrid="full", ): super().__init__( alphas=alphas, @@ -215,4 +238,5 @@ def __init__( verbose=verbose, n_jobs=n_jobs, method=method, + xygrid=xygrid, ) diff --git a/mrinversion/linear_model/utils.py b/mrinversion/linear_model/utils.py new file mode 100644 index 0000000..0ad5d48 --- /dev/null +++ b/mrinversion/linear_model/utils.py @@ -0,0 +1,225 @@ +import numpy as np +from scipy.optimize import minimize + +from mrinversion.linear_model import SmoothLassoCV + + +class CVMinimizer: + def __init__( + self, + inverse_dimension, + compressed_K, + compressed_s, + guess_neglogalpha, + guess_negloglambda, + sigma=0.0, + folds=10, + xygrid="full", + ): + self.inverse_dimension = inverse_dimension + self.compressed_K = compressed_K + self.compressed_s = compressed_s + self.guess = np.array([guess_neglogalpha, guess_negloglambda]) + self.sigma = sigma + self.folds = folds + self.xygrid = xygrid + + @staticmethod + def _cv_min_withprint(loghyperparameters, self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid, + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + print(f"alpha: {-np.log10(alpha)}, lambda: {-np.log10(lam)}, cv: {cv}") + return cv + + @staticmethod + def _cv_min_noprint(loghyperparameters, self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid, + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + return cv + + @staticmethod + def _print_result(this_result, print_steps, prev_result=None): + if prev_result: + if this_result.fun < prev_result.fun: + if this_result.success: + print( + "This minimization improved on the previous result and completed the minimization in the iterations given. " + ) + else: + print( + "This minimization improved on the previous result, but did not finish in the given number of iterations." + ) + else: + if this_result.success: + print("This minimization minimized to a worse solution.") + else: + print( + "This minimization did not improve on the previous result, and did not finish in the given number of iterations." + ) + if print_steps: + print(f"This minimization results:\n\n {this_result}") + else: + if this_result.success: + print("This minimization completed in the iterations given. ") + else: + print( + "This minimization did not finish in the given number of iterations." + ) + if print_steps: + print(f"This minimization results:\n\n {this_result}") + + print("--------------------------------------------------\n") + + def _minimize( + self, guess, tol, maxiter, print_funcevals, print_steps, prev_result=None + ): + if print_funcevals: + simplex_res = minimize( + fun=self._cv_min_withprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={"maxiter": maxiter}, + method="Nelder-Mead", + ) + else: + simplex_res = minimize( + fun=self._cv_min_noprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={"maxiter": maxiter}, + method="Nelder-Mead", + ) + self._print_result( + this_result=simplex_res, print_steps=print_steps, prev_result=prev_result + ) + return simplex_res + + def mincv_close( + self, tol=0.0002, maxiter=500, print_steps=True, print_funcevals=False + ): + simplex_res1 = self._minimize( + guess=self.guess, + tol=tol * 1.5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + ) + + if simplex_res1.success: + simplex_res2 = self._minimize( + guess=simplex_res1.x, + tol=tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + prev_result=simplex_res1, + ) + if simplex_res2.fun < simplex_res1.fun: + return simplex_res2 + else: + return simplex_res1 + else: + return simplex_res1 + + def mincv( + self, + tol=0.0002, + maxiter=100, + print_steps=True, + print_funcevals=False, + guess_bounds=[[2, 10], [2, 10]], + guesses="all", + ): + '''options for guesses: "all", "cartesian", "diagonals"''' + lowval_alpha = np.average([self.guess[0], guess_bounds[0][0]]) + highval_alpha = np.average([self.guess[0], guess_bounds[0][1]]) + lowval_lambda = np.average([self.guess[1], guess_bounds[1][0]]) + highval_lambda = np.average([self.guess[1], guess_bounds[1][1]]) + if guesses == "all": + trythese = np.zeros((9, 2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, highval_lambda] + trythese[2] = [highval_alpha, self.guess[1]] + trythese[3] = [highval_alpha, lowval_lambda] + trythese[4] = [self.guess[0], lowval_lambda] + trythese[5] = [lowval_alpha, lowval_lambda] + trythese[6] = [lowval_alpha, self.guess[1]] + trythese[7] = [lowval_alpha, highval_lambda] + trythese[8] = self.guess + elif guesses == "cartesian": + trythese = np.zeros((5, 2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, self.guess[1]] + trythese[2] = [self.guess[0], lowval_lambda] + trythese[3] = [lowval_alpha, self.guess[1]] + trythese[4] = self.guess + elif guesses == "diagonal": + trythese = np.zeros((5, 2)) + trythese[0] = [highval_alpha, highval_lambda] + trythese[1] = [highval_alpha, lowval_lambda] + trythese[2] = [lowval_alpha, lowval_lambda] + trythese[3] = [lowval_alpha, highval_lambda] + trythese[4] = self.guess + else: + print('choose a valid choice for "guesses"') + return None + print(f"Starting points to try: {trythese}") + first_step = [ + self._minimize( + guess=thisguess, + tol=tol * 5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + ) + for thisguess in trythese + ] + print("-----------------------------------------------------------------") + mins = np.asarray([step.fun for step in first_step]) + min_idx = np.where(mins == mins.min())[0][0] + + final_step = self._minimize( + guess=first_step[min_idx].x, + tol=tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=False, + ) + + self._print_result( + final_step, print_steps=True, prev_result=first_step[min_idx] + ) + + if final_step.fun > first_step[min_idx].fun: + print("Second minimization did not improve.") + return first_step[min_idx] + else: + return final_step From b29a5e50d9d2bb04d5e9d09513704478b6ec06b0 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:42:56 -0400 Subject: [PATCH 09/20] linting --- mrinversion/kernel/base.py | 3 +-- mrinversion/kernel/csa_aniso.py | 5 ---- mrinversion/kernel/nmr.py | 5 +--- mrinversion/kernel/quad_aniso.py | 24 +++++++++++++++++--- mrinversion/linear_model/_base_l1l2.py | 1 - mrinversion/linear_model/linear_inversion.py | 4 ---- mrinversion/linear_model/smooth_lasso.py | 2 -- 7 files changed, 23 insertions(+), 21 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 64d85e4..bae885b 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -2,7 +2,6 @@ import csdmpy as cp import numpy as np -import numpy.ma as ma from mrsimulator.method.lib import BlochDecaySpectrum from .utils import _x_y_to_cq_eta_distribution @@ -165,7 +164,7 @@ def __init__( if number_of_sidebands is None: self.number_of_sidebands = dim.count - def _get_zeta_eta(self, supersampling, eta_bound=1,calc_pos=False): + def _get_zeta_eta(self, supersampling, eta_bound=1, calc_pos=False): """Return zeta and eta coordinates over x-y grid""" if eta_bound == 1 and not calc_pos: zeta, eta = _x_y_to_zeta_eta_distribution( diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index 18b3989..4d47ebd 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -1,13 +1,8 @@ from copy import deepcopy -import csdmpy as cp -import numpy as np from mrsimulator import Simulator from mrsimulator import SpinSystem -from mrsimulator.method import Method -from mrsimulator.method import SpectralEvent from mrsimulator.method.lib import BlochDecaySpectrum -from mrsimulator.utils import get_spectral_dimensions from mrinversion.kernel.base import LineShape diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index 3a1d9e2..c589e70 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -5,9 +5,6 @@ from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA -from mrinversion.kernel.quad_aniso import DAS # NOQA -from mrinversion.kernel.quad_aniso import MQMAS # NOQA -from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA -from mrinversion.kernel.quad_aniso import SL_MQMASnodist + # from mrinversion.kernel.quad_aniso import MQMAS # NOQA diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 3c8cf82..ed22c30 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -73,6 +73,13 @@ def kernel(self, supersampling, eta_bound=1): amp = sim.methods[0].simulation.real + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we " \ + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ + "sign of Cq from other means, you can restrict the grid using " \ + "xygrid='positive' or xygrid='negative." + ) + return self._averaged_kernel(amp, supersampling) @@ -131,7 +138,10 @@ def kernel(self, supersampling, eta_bound=1): amp = sim.methods[0].simulation.real warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + "This kernel is intended to be used with xygrid='mirrored', since we " \ + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ + "sign of Cq from other means, you can restrict the grid using " \ + "xygrid='positive' or xygrid='negative." ) return self._averaged_kernel(amp, supersampling) @@ -211,8 +221,12 @@ def kernel(self, supersampling, eta_bound=1): ) warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + "This kernel is intended to be used with xygrid='mirrored', since we " \ + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ + "sign of Cq from other means, you can restrict the grid using " \ + "xygrid='positive' or xygrid='negative." ) + return self._averaged_kernel(amp, supersampling) @@ -293,8 +307,12 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True): ) warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we cannot distinguish +Cq from -Cq from an NMR experiment. If you know the sign of Cq from other means, you can restrict the grid using xygrid='positive' or xygrid='negative." + "This kernel is intended to be used with xygrid='mirrored', since we " \ + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ + "sign of Cq from other means, you can restrict the grid using " \ + "xygrid='positive' or xygrid='negative." ) + return self._averaged_kernel(amp, supersampling) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 2ac27ea..5b47ba6 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -352,7 +352,6 @@ def __init__( self.f_shape = tuple(item.count for item in inverse_dimension)[::-1] self.xygrid = xygrid - def fit(self, K, s, cv_map_as_csdm=True): def fit(self, K, s, cv_map_as_csdm=True): r"""Fit the model using the coordinate descent method from scikit-learn for all alpha anf lambda values using the `n`-folds cross-validation technique. diff --git a/mrinversion/linear_model/linear_inversion.py b/mrinversion/linear_model/linear_inversion.py index 9900031..25b3d92 100644 --- a/mrinversion/linear_model/linear_inversion.py +++ b/mrinversion/linear_model/linear_inversion.py @@ -25,10 +25,6 @@ def find_optimum_singular_value(s): def TSVD(K): - # print(f'K: {K}') - # if ma.is_masked(K): - # U, S, VT = np.linalg.svd(K.filled(0), full_matrices=False) - # else: U, S, VT = np.linalg.svd(K, full_matrices=False) r = find_optimum_singular_value(S) return U, S, VT, r diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index cd3236d..3ee4b9d 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -1,5 +1,3 @@ -import warnings - from ._base_l1l2 import GeneralL2Lasso from ._base_l1l2 import GeneralL2LassoCV From 38e17d0ac721065af8835d8b55fe570cf22975a3 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:51:06 -0400 Subject: [PATCH 10/20] more linting, remove unneeded if block --- mrinversion/kernel/quad_aniso.py | 26 ++++---- mrinversion/linear_model/_base_l1l2.py | 75 +++++++++--------------- mrinversion/linear_model/smooth_lasso.py | 2 - 3 files changed, 40 insertions(+), 63 deletions(-) diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index ed22c30..3ac6168 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -74,9 +74,9 @@ def kernel(self, supersampling, eta_bound=1): amp = sim.methods[0].simulation.real warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we " \ - "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ - "sign of Cq from other means, you can restrict the grid using " \ + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " "xygrid='positive' or xygrid='negative." ) @@ -138,10 +138,10 @@ def kernel(self, supersampling, eta_bound=1): amp = sim.methods[0].simulation.real warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we " \ - "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ - "sign of Cq from other means, you can restrict the grid using " \ - "xygrid='positive' or xygrid='negative." + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " + "xygrid='positive' or xygrid='negative'." ) return self._averaged_kernel(amp, supersampling) @@ -221,9 +221,9 @@ def kernel(self, supersampling, eta_bound=1): ) warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we " \ - "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ - "sign of Cq from other means, you can restrict the grid using " \ + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " "xygrid='positive' or xygrid='negative." ) @@ -307,9 +307,9 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True): ) warnings.warn( - "This kernel is intended to be used with xygrid='mirrored', since we " \ - "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" \ - "sign of Cq from other means, you can restrict the grid using " \ + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " "xygrid='positive' or xygrid='negative." ) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 5b47ba6..b027c42 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -64,13 +64,11 @@ def __init__( ): if ( - xygrid != "full" - and xygrid != "mirrored" - and xygrid != "positive" - and xygrid != "negative" + xygrid != "full" and xygrid != "mirrored" and xygrid != "positive" and xygrid != "negative" ): raise Exception( - "Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'." + "Please choose a valid option for 'xygrid'. Valid options are 'full', " + "'mirrored', 'positive', and 'negative'." ) self.hyperparameters = {"lambda": lambda1, "alpha": alpha} @@ -208,7 +206,6 @@ def _sol_to_csdm(self, s): if self.xygrid != "full": if self.xygrid == "mirrored": - # print('made it to mirrored') self.f = self.mirror_inversion() elif self.xygrid == "negative": self.f = self.flip_inversion() @@ -322,13 +319,11 @@ def __init__( xygrid="full", ): if ( - xygrid != "full" - and xygrid != "mirrored" - and xygrid != "positive" - and xygrid != "negative" + xygrid != "full" and xygrid != "mirrored" and xygrid != "positive" and xygrid != "negative" ): raise Exception( - "Please choose a valid option for 'xygrid'. Valid options are 'full', 'mirrored', 'positive', and 'negative'." + "Please choose a valid option for 'xygrid'. Valid options are 'full', " + "'mirrored', 'positive', and 'negative'." ) self.cv_alphas = np.asarray(alphas).ravel() @@ -442,46 +437,30 @@ def fit(self, K, s, cv_map_as_csdm=True): # Calculate the solution using the complete data at the optimized lambda and # alpha values - if self.xygrid == "full": - self.opt = GeneralL2Lasso( - alpha=self.hyperparameters["alpha"], - lambda1=self.hyperparameters["lambda"], - max_iterations=self.max_iterations, - tolerance=self.tolerance, - positive=self.positive, - regularizer=self.regularizer, - inverse_dimension=self.inverse_dimension, - method=self.method, - ) - self.opt.fit(K, s) - self.f = self.opt.f - if cv_map_as_csdm: - self.cv_map_to_csdm() - - else: - self.opt = GeneralL2Lasso( - alpha=self.hyperparameters["alpha"], - lambda1=self.hyperparameters["lambda"], - max_iterations=self.max_iterations, - tolerance=self.tolerance, - positive=self.positive, - regularizer=self.regularizer, - inverse_dimension=self.inverse_dimension, - method=self.method, - xygrid=self.xygrid, - ) + + self.opt = GeneralL2Lasso( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + xygrid=self.xygrid, + ) - self.opt.fit(K, s) - self.f = self.opt.f + self.opt.fit(K, s) + self.f = self.opt.f - if self.xygrid == "mirrored": - # print('made it to mirrored') - self.f = self.mirror_inversion() - elif self.xygrid == "negative": - self.f = self.flip_inversion() + if self.xygrid == "mirrored": + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() - if cv_map_as_csdm: - self.cv_map_to_csdm() + if cv_map_as_csdm: + self.cv_map_to_csdm() def cv_map_to_csdm(self): # convert cv_map to csdm diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index 3ee4b9d..297ff72 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -1,8 +1,6 @@ from ._base_l1l2 import GeneralL2Lasso from ._base_l1l2 import GeneralL2LassoCV -import warnings - __author__ = "Deepansh J. Srivastava" __email__ = "srivastava.89@osu.edu" From 4c040c67291c3279285a6fd1a4d23c61124b4218 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:55:10 -0400 Subject: [PATCH 11/20] even more linting --- mrinversion/linear_model/_base_l1l2.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index b027c42..7c8b99e 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -64,10 +64,13 @@ def __init__( ): if ( - xygrid != "full" and xygrid != "mirrored" and xygrid != "positive" and xygrid != "negative" + xygrid != "full" and + xygrid != "mirrored" and + xygrid != "positive" and + xygrid != "negative" ): raise Exception( - "Please choose a valid option for 'xygrid'. Valid options are 'full', " + "Please choose a valid option for 'xygrid'. Valid options are 'full', " "'mirrored', 'positive', and 'negative'." ) @@ -319,7 +322,10 @@ def __init__( xygrid="full", ): if ( - xygrid != "full" and xygrid != "mirrored" and xygrid != "positive" and xygrid != "negative" + xygrid != "full" and + xygrid != "mirrored" and + xygrid != "positive" and + xygrid != "negative" ): raise Exception( "Please choose a valid option for 'xygrid'. Valid options are 'full', " From 544002ab2d0eee3906fad4247101bc19654e1606 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 15:59:07 -0400 Subject: [PATCH 12/20] fixing if statement again, fix more whitespace stuff --- mrinversion/kernel/quad_aniso.py | 1 - mrinversion/linear_model/_base_l1l2.py | 17 ++++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 3ac6168..7c00e57 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -312,7 +312,6 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True): "sign of Cq from other means, you can restrict the grid using " "xygrid='positive' or xygrid='negative." ) - return self._averaged_kernel(amp, supersampling) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 7c8b99e..0e81564 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -64,10 +64,10 @@ def __init__( ): if ( - xygrid != "full" and - xygrid != "mirrored" and - xygrid != "positive" and - xygrid != "negative" + xygrid != "full" + and xygrid != "mirrored" + and xygrid != "positive" + and xygrid != "negative" ): raise Exception( "Please choose a valid option for 'xygrid'. Valid options are 'full', " @@ -322,10 +322,10 @@ def __init__( xygrid="full", ): if ( - xygrid != "full" and - xygrid != "mirrored" and - xygrid != "positive" and - xygrid != "negative" + xygrid != "full" + and xygrid != "mirrored" + and xygrid != "positive" + and xygrid != "negative" ): raise Exception( "Please choose a valid option for 'xygrid'. Valid options are 'full', " @@ -443,7 +443,6 @@ def fit(self, K, s, cv_map_as_csdm=True): # Calculate the solution using the complete data at the optimized lambda and # alpha values - self.opt = GeneralL2Lasso( alpha=self.hyperparameters["alpha"], lambda1=self.hyperparameters["lambda"], From 9b8a2973e8e9a3eb59dd8904b835c76e452a27b6 Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 16:04:21 -0400 Subject: [PATCH 13/20] remove import line --- mrinversion/kernel/base.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index bae885b..b0d7ad9 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -4,7 +4,6 @@ import numpy as np from mrsimulator.method.lib import BlochDecaySpectrum -from .utils import _x_y_to_cq_eta_distribution from .utils import _x_y_to_zeta_eta_distribution __dimension_list__ = (cp.Dimension, cp.LinearDimension, cp.MonotonicDimension) From 6e250f064e3b03695d67658e467b8f1370e3cd2e Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 29 Aug 2025 16:09:03 -0400 Subject: [PATCH 14/20] get rid of more old grid stuff --- mrinversion/kernel/base.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index b0d7ad9..a4f6657 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -176,14 +176,6 @@ def _get_zeta_eta(self, supersampling, eta_bound=1, calc_pos=False): ) return zeta, eta, abundances - def _get_cq_eta(self, supersampling): - """Return zeta and eta coordinates over x-y grid""" - - cq, eta = _x_y_to_cq_eta_distribution( - self.inverse_kernel_dimension, supersampling - ) - return cq, eta - def _check_csdm_dimension(dimensions, dimension_id): if not isinstance(dimensions, (list, *__dimension_list__)): From be620308803da2dc4d199f73a50dedcd7f8f17e8 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 1 Sep 2025 12:16:32 -0400 Subject: [PATCH 15/20] lint --- mrinversion/kernel/nmr.py | 8 ++++---- mrinversion/kernel/utils.py | 3 +-- mrinversion/linear_model/linear_inversion.py | 1 + mrinversion/linear_model/utils.py | 9 ++++++--- mrinversion/utils.py | 19 +++++++++---------- 5 files changed, 21 insertions(+), 19 deletions(-) diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index c589e70..c648566 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,10 +1,10 @@ -from mrinversion.kernel.quad_aniso import DAS # NOQA -from mrinversion.kernel.quad_aniso import MQMAS # NOQA -from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA -from mrinversion.kernel.quad_aniso import SL_MQMASnodist from mrinversion.kernel.csa_aniso import MAF # NOQA from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA +from mrinversion.kernel.quad_aniso import DAS # NOQA +from mrinversion.kernel.quad_aniso import MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMASnodist # NOQA # from mrinversion.kernel.quad_aniso import MQMAS # NOQA diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index 5a8c792..b882758 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -132,7 +132,7 @@ def zeta_eta_to_x_y(zeta, eta): return x.ravel(), y.ravel() -def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound, calc_pos=False): +def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound=1, calc_pos=False): """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) @@ -150,7 +150,6 @@ def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound, calc_pos=False ) return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound, calc_pos) - return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound, calc_pos) # def _x_y_to_cq_eta_distribution(grid, supersampling): diff --git a/mrinversion/linear_model/linear_inversion.py b/mrinversion/linear_model/linear_inversion.py index 25b3d92..374de5c 100644 --- a/mrinversion/linear_model/linear_inversion.py +++ b/mrinversion/linear_model/linear_inversion.py @@ -1,4 +1,5 @@ import numpy as np + # import numpy.ma as ma __author__ = "Deepansh J. Srivastava" diff --git a/mrinversion/linear_model/utils.py b/mrinversion/linear_model/utils.py index c5bc630..b8bf1c7 100644 --- a/mrinversion/linear_model/utils.py +++ b/mrinversion/linear_model/utils.py @@ -67,18 +67,21 @@ def _print_result(this_result, print_steps, prev_result=None): if this_result.fun < prev_result.fun: if this_result.success: print( - "This minimization improved on the previous result and completed the minimization in the iterations given. " + "This minimization improved on the previous result and \n" + "completed the minimization in the iterations given. " ) else: print( - "This minimization improved on the previous result, but did not finish in the given number of iterations." + "This minimization improved on the previous result, but did \n" + "not finish in the given number of iterations." ) else: if this_result.success: print("This minimization minimized to a worse solution.") else: print( - "This minimization did not improve on the previous result, and did not finish in the given number of iterations." + "This minimization did not improve on the previous result, \n" + "and did not finish in the given number of iterations." ) if print_steps: print(f"This minimization results:\n\n {this_result}") diff --git a/mrinversion/utils.py b/mrinversion/utils.py index c8ce5d6..1a9e97c 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -36,7 +36,7 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): if len(csdm_object.x) > 2: extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) data.shape = (extra_dims, data.shape[-2], data.shape[-1]) - + reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) dx = reg_x[1] - reg_x[0] dy = reg_y[1] - reg_y[0] @@ -142,10 +142,10 @@ def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): avg_range_eta = (np.arange(n) - (n - 1) / 2) * d_eta / n # print(avg_range_cq) # print(avg_range_eta) - for i,cq_item in enumerate(avg_range_cq): - for j,eta_item in enumerate(avg_range_eta): + for i, cq_item in enumerate(avg_range_cq): + for j, eta_item in enumerate(avg_range_eta): # print(i,j) - cq__ = (reg_cq + cq_item) + cq__ = reg_cq + cq_item # print(np.abs(cq__).shape) eta__ = np.abs(reg_eta + eta_item) cq_, eta_ = np.meshgrid(cq__, eta__) @@ -158,21 +158,20 @@ def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): # print(eta_) # print() - theta = np.ones(eta_.shape)*1e10 + theta = np.ones(eta_.shape) * 1e10 # print(theta.shape) index = np.where(cq_ > 0) # print(f'pos index: {index}') - theta[index] = np.pi/2 * (1-eta_[index] / 2.0) + theta[index] = np.pi / 2 * (1 - eta_[index] / 2.0) index = np.where(cq_ <= 0) # print(f'neg index: {index}') - - theta[index] = np.pi/4.0 * eta_[index] + + theta[index] = np.pi / 4.0 * eta_[index] # print(cq_) # print(np.where(theta * 180/np.pi==1e10)) # print(list(zip(cq_, theta * 180/np.pi))) - xq_grid = np.abs(cq_) * np.sin(theta) yq_grid = np.abs(cq_) * np.cos(theta) @@ -621,4 +620,4 @@ def plot_3d_( ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) \ No newline at end of file + ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) From 401780fff4d6ad7d96f858ee2e1d322ccd6729b3 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 1 Sep 2025 12:33:07 -0400 Subject: [PATCH 16/20] lint --- .flake8 | 2 +- docs/conf.py | 2 +- mrinversion/kernel/utils.py | 2 +- mrinversion/linear_model/_base_l1l2.py | 20 +++++--------------- mrinversion/linear_model/utils.py | 8 +++----- 5 files changed, 11 insertions(+), 23 deletions(-) diff --git a/.flake8 b/.flake8 index 0909542..7b04ae9 100644 --- a/.flake8 +++ b/.flake8 @@ -1,6 +1,6 @@ [flake8] ignore = E203 -exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py +exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py, examples/* filename = *.pyx, *.px*, *py max-line-length = 88 max-complexity = 10 diff --git a/docs/conf.py b/docs/conf.py index be2e5f0..5e49602 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -248,7 +248,7 @@ # Theme options html_logo = "_static/mrinversion.png" html_style = "style.css" -html_title = f"mrinversion:doc v{__version__}" +html_title = f"mrinversion: doc v{__version__}" html_last_updated_fmt = "" # html_logo = "mrinversion" html_sidebars = { diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index b882758..fea9ace 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -34,7 +34,7 @@ def x_y_to_zeta_eta(x, y): y = y.value if x_unit != y_unit: raise ValueError( - f"x and y must have same dimensionality; x ({x_unit}) != y ({y_unit})." + f"x and y must have same dimensionality: x ({x_unit}) != y ({y_unit})." ) zeta = np.sqrt(x**2 + y**2) # + offset diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 0e81564..d7ee368 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -14,6 +14,7 @@ __author__ = "Deepansh J. Srivastava" __email__ = "srivastava.89@osu.edu" +ALLOWED_GRIDS = ("full", "mirrored", "positive", "negative") class GeneralL2Lasso: @@ -62,14 +63,8 @@ def __init__( method="gradient_decent", xygrid="full", ): - - if ( - xygrid != "full" - and xygrid != "mirrored" - and xygrid != "positive" - and xygrid != "negative" - ): - raise Exception( + if xygrid not in ALLOWED_GRIDS: + raise KeyError( "Please choose a valid option for 'xygrid'. Valid options are 'full', " "'mirrored', 'positive', and 'negative'." ) @@ -321,13 +316,8 @@ def __init__( method="gradient_decent", xygrid="full", ): - if ( - xygrid != "full" - and xygrid != "mirrored" - and xygrid != "positive" - and xygrid != "negative" - ): - raise Exception( + if xygrid not in ALLOWED_GRIDS: + raise KeyError( "Please choose a valid option for 'xygrid'. Valid options are 'full', " "'mirrored', 'positive', and 'negative'." ) diff --git a/mrinversion/linear_model/utils.py b/mrinversion/linear_model/utils.py index b8bf1c7..a65d071 100644 --- a/mrinversion/linear_model/utils.py +++ b/mrinversion/linear_model/utils.py @@ -84,16 +84,14 @@ def _print_result(this_result, print_steps, prev_result=None): "and did not finish in the given number of iterations." ) if print_steps: - print(f"This minimization results:\n\n {this_result}") + print(f"This minimization results: \n\n{this_result}") else: if this_result.success: print("This minimization completed in the iterations given. ") else: - print( - "This minimization did not finish in the given number of iterations." - ) + print("This minimization did not finish in given number of iterations.") if print_steps: - print(f"This minimization results:\n\n {this_result}") + print(f"This minimization results: \n\n{this_result}") print("--------------------------------------------------\n") From c4e8409567421c56821aa79a297102713c833431 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 1 Sep 2025 12:34:57 -0400 Subject: [PATCH 17/20] lint --- mrinversion/linear_model/_base_l1l2.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index d7ee368..65c8932 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -202,11 +202,10 @@ def _sol_to_csdm(self, s): self.f.dimensions[1] = s.dimensions[1] self.f.dimensions[0] = self.inverse_dimension[0] - if self.xygrid != "full": - if self.xygrid == "mirrored": - self.f = self.mirror_inversion() - elif self.xygrid == "negative": - self.f = self.flip_inversion() + if self.xygrid == "mirrored": + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() def flip_inversion(self): csdm_obj = self.f.copy() From d798d0a8796563c704be7b762709d6486476c37e Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 7 Nov 2025 13:12:14 -0500 Subject: [PATCH 18/20] add nQ, add comment line to avoid merge conflict --- mrinversion/kernel/base.py | 1 + mrinversion/kernel/quad_aniso.py | 36 +++++++++++--------------------- 2 files changed, 13 insertions(+), 24 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index a4f6657..ad16f9c 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -1,3 +1,4 @@ +"""Base classes for kernel generation.""" from copy import deepcopy import csdmpy as cp diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 7c00e57..7f9a8eb 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -171,7 +171,7 @@ def __init__( self.exp_dict = exp_dict self.anisotropic_dimension = anisotropic_dimension - def kernel(self, supersampling, eta_bound=1): + def kernel(self, supersampling, eta_bound=1,nQ=3): import sys sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") @@ -254,7 +254,7 @@ def __init__( self.exp_dict = exp_dict self.anisotropic_dimension = anisotropic_dimension - def kernel(self, supersampling, eta_bound=1, cq_posneg=True): + def kernel(self, supersampling, eta_bound=1, cq_posneg=True,n_quantum=3): import sys sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") @@ -265,32 +265,19 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True): isotope = self.method_args["channels"][0] - if eta_bound == 1 and cq_posneg: - Cq, eta = self._get_zeta_eta(supersampling, eta_bound) - else: - Cq, eta, abundances = self._get_zeta_eta( - supersampling, eta_bound, cq_posneg - ) + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound,calc_pos=True) - if eta_bound == 1: - spin_systems = [ - SpinSystem( - sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] - ) - for cq_, e in zip(Cq, eta) - ] - else: - spin_systems = [ - SpinSystem( - sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], - abundance=abun, - ) - for cq_, e, abun in zip(Cq, eta, abundances) - ] + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) spec_dim = get_spectral_dimensions(obj) - + operators=smsim.build_spin_matrices(nucleus=isotope, n_quantum=n_quantum) amp = np.asarray( [ smsim.simulate_onesite_lineshape( @@ -301,6 +288,7 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True): contribs="c0_c4", return_array=True, distorted=True, + operators=operators ) for mysys in spin_systems ] From 9dd7c16a676566f179c7cbc89643717d909e746f Mon Sep 17 00:00:00 2001 From: mccarthy677 Date: Fri, 7 Nov 2025 13:13:28 -0500 Subject: [PATCH 19/20] finishing merge? --- .flake8 | 18 +- .github/workflows/codeql-analysis.yml | 116 +- CHANGELOG | 60 +- docs/conf.py | 724 +++++----- mrinversion/__init__.py | 2 +- mrinversion/kernel/__init__.py | 6 +- mrinversion/linear_model/linear_inversion.py | 92 +- mrinversion/tests/fista_intergration_test.py | 110 +- mrinversion/utils.py | 1246 +++++++++--------- setup.cfg | 38 +- setup.py | 136 +- 11 files changed, 1274 insertions(+), 1274 deletions(-) diff --git a/.flake8 b/.flake8 index 7b04ae9..1b6de64 100644 --- a/.flake8 +++ b/.flake8 @@ -1,9 +1,9 @@ -[flake8] -ignore = E203 -exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py, examples/* -filename = *.pyx, *.px*, *py -max-line-length = 88 -max-complexity = 10 -select = C,E,F,W,N8 -count = True -statistics = True +[flake8] +ignore = E203 +exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py, examples/* +filename = *.pyx, *.px*, *py +max-line-length = 88 +max-complexity = 10 +select = C,E,F,W,N8 +count = True +statistics = True diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index 84cb6c8..9cf31a8 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -1,58 +1,58 @@ -# For most projects, this workflow file will not need changing; you simply need -# to commit it to your repository. -# -# You may wish to alter this file to override the set of languages analyzed, -# or to provide custom queries or build logic. -name: "CodeQL" - -on: - # At the end of every day - schedule: - - cron: '0 0 * * *' - -jobs: - analyze: - name: Analyze - runs-on: ubuntu-latest - - strategy: - fail-fast: false - matrix: - # Override automatic language detection by changing the below list - # Supported options are ['csharp', 'cpp', 'go', 'java', 'javascript', 'python'] - language: ['python'] - # Learn more... - # https://docs.github.com/en/github/finding-security-vulnerabilities-and-errors-in-your-code/configuring-code-scanning#overriding-automatic-language-detection - - steps: - - name: Checkout repository - uses: actions/checkout@v5 - - # Initializes the CodeQL tools for scanning. - - name: Initialize CodeQL - uses: github/codeql-action/init@v3 - with: - languages: ${{ matrix.language }} - # If you wish to specify custom queries, you can do so here or in a config file. - # By default, queries listed here will override any specified in a config file. - # Prefix the list here with "+" to use these queries and those in the config file. - # queries: ./path/to/local/query, your-org/your-repo/queries@main - - # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). - # If this step fails, then you should remove it and run the build manually (see below) - - name: Autobuild - uses: github/codeql-action/autobuild@v3 - - # ℹ️ Command-line programs to run using the OS shell. - # 📚 https://git.io/JvXDl - - # ✏️ If the Autobuild fails above, remove it and uncomment the following three lines - # and modify them (or add more) to build your code if your project - # uses a compiled language - - #- run: | - # make bootstrap - # make release - - - name: Perform CodeQL Analysis - uses: github/codeql-action/analyze@v3 +# For most projects, this workflow file will not need changing; you simply need +# to commit it to your repository. +# +# You may wish to alter this file to override the set of languages analyzed, +# or to provide custom queries or build logic. +name: "CodeQL" + +on: + # At the end of every day + schedule: + - cron: '0 0 * * *' + +jobs: + analyze: + name: Analyze + runs-on: ubuntu-latest + + strategy: + fail-fast: false + matrix: + # Override automatic language detection by changing the below list + # Supported options are ['csharp', 'cpp', 'go', 'java', 'javascript', 'python'] + language: ['python'] + # Learn more... + # https://docs.github.com/en/github/finding-security-vulnerabilities-and-errors-in-your-code/configuring-code-scanning#overriding-automatic-language-detection + + steps: + - name: Checkout repository + uses: actions/checkout@v5 + + # Initializes the CodeQL tools for scanning. + - name: Initialize CodeQL + uses: github/codeql-action/init@v3 + with: + languages: ${{ matrix.language }} + # If you wish to specify custom queries, you can do so here or in a config file. + # By default, queries listed here will override any specified in a config file. + # Prefix the list here with "+" to use these queries and those in the config file. + # queries: ./path/to/local/query, your-org/your-repo/queries@main + + # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). + # If this step fails, then you should remove it and run the build manually (see below) + - name: Autobuild + uses: github/codeql-action/autobuild@v3 + + # ℹ️ Command-line programs to run using the OS shell. + # 📚 https://git.io/JvXDl + + # ✏️ If the Autobuild fails above, remove it and uncomment the following three lines + # and modify them (or add more) to build your code if your project + # uses a compiled language + + #- run: | + # make bootstrap + # make release + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v3 diff --git a/CHANGELOG b/CHANGELOG index ea9f43a..d353611 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,30 +1,30 @@ -v0.3.1 ------- - -What's new! -''''''''''' - -- Simplified plot_3d function -- Replace fortran code with numba jit python version. - -v0.3.0 ------- - -What's new! -''''''''''' - -- Added T1, T2 relaxation inversion kernel. - -v0.2.0 ------- - -What's new! -''''''''''' - -- Added :func:`~mrinversion.utils.to_Haeberlen_grid` function to convert the 3D :math:`\rho(\delta_\text{iso}, x, y)` - distribution to :math:`\rho(\delta_\text{iso}, \zeta_\sigma, \eta_\sigma)` distribution. - -Changes -''''''' - -- Update code to comply with updated csdmpy library. +v0.3.1 +------ + +What's new! +''''''''''' + +- Simplified plot_3d function +- Replace fortran code with numba jit python version. + +v0.3.0 +------ + +What's new! +''''''''''' + +- Added T1, T2 relaxation inversion kernel. + +v0.2.0 +------ + +What's new! +''''''''''' + +- Added :func:`~mrinversion.utils.to_Haeberlen_grid` function to convert the 3D :math:`\rho(\delta_\text{iso}, x, y)` + distribution to :math:`\rho(\delta_\text{iso}, \zeta_\sigma, \eta_\sigma)` distribution. + +Changes +''''''' + +- Update code to comply with updated csdmpy library. diff --git a/docs/conf.py b/docs/conf.py index 5e49602..f52a298 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,362 +1,362 @@ -# Configuration file for the Sphinx documentation builder. -# -# This file does only contain a selection of the most common options. For a -# full list see the documentation: -# http://www.sphinx-doc.org/en/master/config -# -- Path setup -------------------------------------------------------------- -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# -import os -import sys -import warnings - -from sphinx_gallery.sorting import ExplicitOrder -from sphinx_gallery.sorting import FileNameSortKey - -sys.path.insert(0, os.path.abspath("../..")) - -# -- Project information ----------------------------------------------------- -project = "mrinversion" -copyright = "2020, Deepansh J. Srivastava" -author = "Deepansh J. Srivastava" - -path = os.path.split(__file__)[0] -# get version number from the file -with open(os.path.join(path, "../mrinversion/__init__.py")) as f: - for line in f.readlines(): - if "__version__" in line: - before_keyword, keyword, after_keyword = line.partition("=") - __version__ = after_keyword.strip()[1:-1] - -# The short X.Y version -version = __version__ -# The full version, including alpha/beta/rc tags -release = __version__ - - -# -- General configuration --------------------------------------------------- - -# If your documentation needs a minimal Sphinx version, state it here. -needs_sphinx = "4.0" - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. -extensions = [ - "sphinx.ext.autodoc", - "matplotlib.sphinxext.plot_directive", - "sphinx.ext.doctest", - "sphinx.ext.mathjax", - "sphinx.ext.autosummary", - "sphinx.ext.viewcode", - "sphinx.ext.napoleon", - "sphinx_copybutton", - "sphinxjp.themes.basicstrap", - "sphinx_gallery.gen_gallery", - "sphinx.ext.intersphinx", - "sphinx_tabs.tabs", -] - -autosummary_generate = True - -# ---------------------------------------------------------------------------- # -# Plot directive config # -# ---------------------------------------------------------------------------- # -plot_html_show_source_link = False -plot_rcparams = { - "font.size": 10, - "font.weight": "light", - "font.family": "sans-serif", - "font.sans-serif": "Helvetica", -} - -# ---------------------------------------------------------------------------- # -# Sphinx Gallery config # -# ---------------------------------------------------------------------------- # - -# filter sphinx matplotlib warning -warnings.filterwarnings( - "ignore", - category=UserWarning, - message="Matplotlib is currently using agg, which is a" - " non-GUI backend, so cannot show the figure.", -) - -# numfig config -numfig = True -numfig_secnum_depth = 1 -numfig_format = {"figure": "Figure %s", "table": "Table %s", "code-block": "Listing %s"} - -# math -math_number_all = True - -# sphinx gallery config -sphinx_gallery_conf = { - "examples_dirs": ["../examples"], # path to example scripts - "remove_config_comments": True, - "gallery_dirs": ["galley_examples"], # path to gallery generated output - "within_subsection_order": FileNameSortKey, - "subsection_order": ExplicitOrder( - [ - "../examples/synthetic", - "../examples/sideband", - "../examples/MAF", - "../examples/relaxation", - ] - ), - "reference_url": { - # The module you locally document uses None - "mrinversion": None, - }, - "first_notebook_cell": ( - "# This cell is added by sphinx-gallery\n\n" - "%matplotlib inline\n\n" - "import mrinversion\n" - "print(f'You are using mrinversion v{mrinversion.__version__}')" - ), - # "binder": { - # # Required keys - # "org": "DeepanshS", - # "repo": "mrinversion", - # "branch": "master", - # "binderhub_url": "https://mybinder.org", - # "dependencies": "../requirements.txt", - # # Optional keys - # "filepath_prefix": "docs/_build/html", - # "notebooks_dir": "../../notebooks", - # "use_jupyter_lab": True, - # }, -} - -intersphinx_mapping = { - "matplotlib": ("https://matplotlib.org", None), - "numpy": ("https://numpy.org/doc/stable/", None), - "csdmpy": ("https://csdmpy.readthedocs.io/en/latest/", None), - "astropy": ("https://docs.astropy.org/en/stable/", None), -} - -copybutton_prompt_text = ">>> |\\$ |\\[\\d*\\]: |\\.\\.\\.: |[.][.][.] " -copybutton_prompt_is_regexp = True - -# ---------------------------------------------------------------------------- # - -# Add any paths that contain templates here, relative to this directory. -templates_path = ["_templates"] - -# The suffix(es) of source filenames. -# You can specify multiple suffix as a list of string: -source_suffix = ".rst" - -# The master toc-tree document. -master_doc = "index" - -# autodoc mock modules -autodoc_mock_imports = [] - -# The language for content autogenerated by Sphinx. Refer to documentation -# for a list of supported languages. -# -# This is also used if you do content translation via gettext catalogs. -# Usually you set "language" from the command line for these cases. -language = "en" - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -# This pattern also affects html_static_path and html_extra_path. -exclude_patterns = ["_build", "**.ipynb_checkpoints", "Thumbs.db", ".DS_Store"] - -# The name of the Pygments (syntax highlighting) style to use. -pygments_style = None - - -# -- Options for HTML output ------------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# - -# Some html_theme options are 'alabaster', 'bootstrap', 'sphinx_rtd_theme', -# 'classic', 'basicstrap' -html_theme = "basicstrap" - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. -# -html_theme_options = { - # Set the lang attribute of the html tag. Defaults to 'en' - "lang": "en", - # Disable showing the sidebar. Defaults to 'false' - "nosidebar": False, - # Show header searchbox. Defaults to false. works only "nosidebar=True", - "header_searchbox": False, - # Put the sidebar on the right side. Defaults to false. - "rightsidebar": False, - # Set the width of the sidebar. Defaults to 3 - "sidebar_span": 3, - # Fix navbar to top of screen. Defaults to true - "nav_fixed_top": True, - # Fix the width of the sidebar. Defaults to false - "nav_fixed": True, - # Set the width of the sidebar. Defaults to '900px' - "nav_width": "300px", - # Fix the width of the content area. Defaults to false - "content_fixed": False, - # Set the width of the content area. Defaults to '900px' - "content_width": "900px", - # Fix the width of the row. Defaults to false - "row_fixed": False, - # Disable the responsive design. Defaults to false - "noresponsive": False, - # Disable the responsive footer relbar. Defaults to false - "noresponsiverelbar": False, - # Disable flat design. Defaults to false. - # Works only "bootstrap_version = 3" - "noflatdesign": False, - # Enable Google Web Font. Defaults to false - "googlewebfont": False, - # Set the URL of Google Web Font's CSS. - # Defaults to 'http://fonts.googleapis.com/css?family=Text+Me+One' - # "googlewebfont_url": "http://fonts.googleapis.com/css?family=Roboto", # NOQA - # Set the Style of Google Web Font's CSS. - # Defaults to "font-family: 'Text Me One', sans-serif;" - # "googlewebfont_style": u"font-family: 'Roboto' Regular;", # font-size: 1.5em", - # Set 'navbar-inverse' attribute to header navbar. Defaults to false. - "header_inverse": True, - # Set 'navbar-inverse' attribute to relbar navbar. Defaults to false. - "relbar_inverse": True, - # Enable inner theme by Bootswatch. Defaults to false - "inner_theme": False, - # Set the name of inner theme. Defaults to 'bootswatch-simplex' - # "inner_theme_name": "bootswatch-simplex", - # Select Twitter bootstrap version 2 or 3. Defaults to '3' - "bootstrap_version": "3", - # Show "theme preview" button in header navbar. Defaults to false. - "theme_preview": False, - # Set the Size of Heading text. Defaults to None - # "h1_size": "3.0em", - # "h2_size": "2.6em", - # "h3_size": "2.2em", - # "h4_size": "1.8em", - # "h5_size": "1.9em", - # "h6_size": "1.1em", -} - - -# Theme options -html_logo = "_static/mrinversion.png" -html_style = "style.css" -html_title = f"mrinversion: doc v{__version__}" -html_last_updated_fmt = "" -# html_logo = "mrinversion" -html_sidebars = { - "**": ["searchbox.html", "globaltoc.html"], - "using/windows": ["searchbox.html", "windowssidebar.html"], -} - - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ["_static"] - - -# -- Options for HTMLHelp output --------------------------------------------- - -# Output file base name for HTML help builder. -htmlhelp_basename = "Mrinversion doc" - -# -- Options for LaTeX output ------------------------------------------------ -latex_engine = "xelatex" -# latex_logo = "_static/csdmpy.png" -latex_show_pagerefs = True - -latex_elements = { - # The paper size ('letterpaper' or 'a4paper'). - # - "papersize": "letterpaper", - # The font size ('10pt', '11pt' or '12pt'). - # - "pointsize": "9pt", - "fontenc": r"\usepackage[utf8]{inputenc}", - "geometry": r"\usepackage[vmargin=2.5cm, hmargin=2cm]{geometry}", - # "fncychap": r"\usepackage[Rejne]{fncychap}", - # Additional stuff for the LaTeX preamble. - "preamble": r""" - \usepackage[T1]{fontenc} - \usepackage{amsfonts, amsmath, amssymb, mathbbol} - \usepackage{graphicx} - \usepackage{setspace} - \singlespacing - - \usepackage{fancyhdr} - \pagestyle{fancy} - \fancyhf{} - \fancyhead[L]{ - \ifthenelse{\isodd{\value{page}}}{ \small \nouppercase{\leftmark} }{} - } - \fancyhead[R]{ - \ifthenelse{\isodd{\value{page}}}{}{ \small \nouppercase{\rightmark} } - } - \fancyfoot[CO, CE]{\thepage} - """, - # Latex figure (float) alignment - # - "figure_align": "htbp", -} - -# Grouping the document tree into LaTeX files. List of tuples -# (source start file, target name, title, -# author, documentclass [howto, manual, or own class]). -latex_documents = [ - (master_doc, "mrinversion.tex", "mrinversion Documentation", author, "manual") -] - - -# -- Options for manual page output ------------------------------------------ - -# One entry per manual page. List of tuples -# (source start file, name, description, authors, manual section). -man_pages = [(master_doc, "mrinversion", "mrinversion Documentation", [author], 1)] - - -# -- Options for Texinfo output ---------------------------------------------- - -# Grouping the document tree into Texinfo files. List of tuples -# (source start file, target name, title, author, -# dir menu entry, description, category) -texinfo_documents = [ - ( - master_doc, - "Mrinversion", - "Mrinversion Documentation", - author, - "Mrinversion", - "Statistical learning of tensor distribution from NMR anisotropic spectra", - "Miscellaneous", - ) -] - - -# -- Options for Epub output ------------------------------------------------- - -# Bibliographic Dublin Core info. -epub_title = project - -# The unique identifier of the text. This can be a ISBN number -# or the project homepage. -# -# epub_identifier = '' - -# A unique identification for the text. -# -# epub_uid = '' - -# A list of files that should not be packed into the epub file. -epub_exclude_files = ["search.html", "_static/style.css"] - - -def setup(app): - app.add_css_file("style.css") +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config +# -- Path setup -------------------------------------------------------------- +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +import warnings + +from sphinx_gallery.sorting import ExplicitOrder +from sphinx_gallery.sorting import FileNameSortKey + +sys.path.insert(0, os.path.abspath("../..")) + +# -- Project information ----------------------------------------------------- +project = "mrinversion" +copyright = "2020, Deepansh J. Srivastava" +author = "Deepansh J. Srivastava" + +path = os.path.split(__file__)[0] +# get version number from the file +with open(os.path.join(path, "../mrinversion/__init__.py")) as f: + for line in f.readlines(): + if "__version__" in line: + before_keyword, keyword, after_keyword = line.partition("=") + __version__ = after_keyword.strip()[1:-1] + +# The short X.Y version +version = __version__ +# The full version, including alpha/beta/rc tags +release = __version__ + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +needs_sphinx = "4.0" + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + "sphinx.ext.autodoc", + "matplotlib.sphinxext.plot_directive", + "sphinx.ext.doctest", + "sphinx.ext.mathjax", + "sphinx.ext.autosummary", + "sphinx.ext.viewcode", + "sphinx.ext.napoleon", + "sphinx_copybutton", + "sphinxjp.themes.basicstrap", + "sphinx_gallery.gen_gallery", + "sphinx.ext.intersphinx", + "sphinx_tabs.tabs", +] + +autosummary_generate = True + +# ---------------------------------------------------------------------------- # +# Plot directive config # +# ---------------------------------------------------------------------------- # +plot_html_show_source_link = False +plot_rcparams = { + "font.size": 10, + "font.weight": "light", + "font.family": "sans-serif", + "font.sans-serif": "Helvetica", +} + +# ---------------------------------------------------------------------------- # +# Sphinx Gallery config # +# ---------------------------------------------------------------------------- # + +# filter sphinx matplotlib warning +warnings.filterwarnings( + "ignore", + category=UserWarning, + message="Matplotlib is currently using agg, which is a" + " non-GUI backend, so cannot show the figure.", +) + +# numfig config +numfig = True +numfig_secnum_depth = 1 +numfig_format = {"figure": "Figure %s", "table": "Table %s", "code-block": "Listing %s"} + +# math +math_number_all = True + +# sphinx gallery config +sphinx_gallery_conf = { + "examples_dirs": ["../examples"], # path to example scripts + "remove_config_comments": True, + "gallery_dirs": ["galley_examples"], # path to gallery generated output + "within_subsection_order": FileNameSortKey, + "subsection_order": ExplicitOrder( + [ + "../examples/synthetic", + "../examples/sideband", + "../examples/MAF", + "../examples/relaxation", + ] + ), + "reference_url": { + # The module you locally document uses None + "mrinversion": None, + }, + "first_notebook_cell": ( + "# This cell is added by sphinx-gallery\n\n" + "%matplotlib inline\n\n" + "import mrinversion\n" + "print(f'You are using mrinversion v{mrinversion.__version__}')" + ), + # "binder": { + # # Required keys + # "org": "DeepanshS", + # "repo": "mrinversion", + # "branch": "master", + # "binderhub_url": "https://mybinder.org", + # "dependencies": "../requirements.txt", + # # Optional keys + # "filepath_prefix": "docs/_build/html", + # "notebooks_dir": "../../notebooks", + # "use_jupyter_lab": True, + # }, +} + +intersphinx_mapping = { + "matplotlib": ("https://matplotlib.org", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "csdmpy": ("https://csdmpy.readthedocs.io/en/latest/", None), + "astropy": ("https://docs.astropy.org/en/stable/", None), +} + +copybutton_prompt_text = ">>> |\\$ |\\[\\d*\\]: |\\.\\.\\.: |[.][.][.] " +copybutton_prompt_is_regexp = True + +# ---------------------------------------------------------------------------- # + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +source_suffix = ".rst" + +# The master toc-tree document. +master_doc = "index" + +# autodoc mock modules +autodoc_mock_imports = [] + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = "en" + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ["_build", "**.ipynb_checkpoints", "Thumbs.db", ".DS_Store"] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = None + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# + +# Some html_theme options are 'alabaster', 'bootstrap', 'sphinx_rtd_theme', +# 'classic', 'basicstrap' +html_theme = "basicstrap" + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +html_theme_options = { + # Set the lang attribute of the html tag. Defaults to 'en' + "lang": "en", + # Disable showing the sidebar. Defaults to 'false' + "nosidebar": False, + # Show header searchbox. Defaults to false. works only "nosidebar=True", + "header_searchbox": False, + # Put the sidebar on the right side. Defaults to false. + "rightsidebar": False, + # Set the width of the sidebar. Defaults to 3 + "sidebar_span": 3, + # Fix navbar to top of screen. Defaults to true + "nav_fixed_top": True, + # Fix the width of the sidebar. Defaults to false + "nav_fixed": True, + # Set the width of the sidebar. Defaults to '900px' + "nav_width": "300px", + # Fix the width of the content area. Defaults to false + "content_fixed": False, + # Set the width of the content area. Defaults to '900px' + "content_width": "900px", + # Fix the width of the row. Defaults to false + "row_fixed": False, + # Disable the responsive design. Defaults to false + "noresponsive": False, + # Disable the responsive footer relbar. Defaults to false + "noresponsiverelbar": False, + # Disable flat design. Defaults to false. + # Works only "bootstrap_version = 3" + "noflatdesign": False, + # Enable Google Web Font. Defaults to false + "googlewebfont": False, + # Set the URL of Google Web Font's CSS. + # Defaults to 'http://fonts.googleapis.com/css?family=Text+Me+One' + # "googlewebfont_url": "http://fonts.googleapis.com/css?family=Roboto", # NOQA + # Set the Style of Google Web Font's CSS. + # Defaults to "font-family: 'Text Me One', sans-serif;" + # "googlewebfont_style": u"font-family: 'Roboto' Regular;", # font-size: 1.5em", + # Set 'navbar-inverse' attribute to header navbar. Defaults to false. + "header_inverse": True, + # Set 'navbar-inverse' attribute to relbar navbar. Defaults to false. + "relbar_inverse": True, + # Enable inner theme by Bootswatch. Defaults to false + "inner_theme": False, + # Set the name of inner theme. Defaults to 'bootswatch-simplex' + # "inner_theme_name": "bootswatch-simplex", + # Select Twitter bootstrap version 2 or 3. Defaults to '3' + "bootstrap_version": "3", + # Show "theme preview" button in header navbar. Defaults to false. + "theme_preview": False, + # Set the Size of Heading text. Defaults to None + # "h1_size": "3.0em", + # "h2_size": "2.6em", + # "h3_size": "2.2em", + # "h4_size": "1.8em", + # "h5_size": "1.9em", + # "h6_size": "1.1em", +} + + +# Theme options +html_logo = "_static/mrinversion.png" +html_style = "style.css" +html_title = f"mrinversion: doc v{__version__}" +html_last_updated_fmt = "" +# html_logo = "mrinversion" +html_sidebars = { + "**": ["searchbox.html", "globaltoc.html"], + "using/windows": ["searchbox.html", "windowssidebar.html"], +} + + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = "Mrinversion doc" + +# -- Options for LaTeX output ------------------------------------------------ +latex_engine = "xelatex" +# latex_logo = "_static/csdmpy.png" +latex_show_pagerefs = True + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + "papersize": "letterpaper", + # The font size ('10pt', '11pt' or '12pt'). + # + "pointsize": "9pt", + "fontenc": r"\usepackage[utf8]{inputenc}", + "geometry": r"\usepackage[vmargin=2.5cm, hmargin=2cm]{geometry}", + # "fncychap": r"\usepackage[Rejne]{fncychap}", + # Additional stuff for the LaTeX preamble. + "preamble": r""" + \usepackage[T1]{fontenc} + \usepackage{amsfonts, amsmath, amssymb, mathbbol} + \usepackage{graphicx} + \usepackage{setspace} + \singlespacing + + \usepackage{fancyhdr} + \pagestyle{fancy} + \fancyhf{} + \fancyhead[L]{ + \ifthenelse{\isodd{\value{page}}}{ \small \nouppercase{\leftmark} }{} + } + \fancyhead[R]{ + \ifthenelse{\isodd{\value{page}}}{}{ \small \nouppercase{\rightmark} } + } + \fancyfoot[CO, CE]{\thepage} + """, + # Latex figure (float) alignment + # + "figure_align": "htbp", +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, "mrinversion.tex", "mrinversion Documentation", author, "manual") +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [(master_doc, "mrinversion", "mrinversion Documentation", [author], 1)] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + ( + master_doc, + "Mrinversion", + "Mrinversion Documentation", + author, + "Mrinversion", + "Statistical learning of tensor distribution from NMR anisotropic spectra", + "Miscellaneous", + ) +] + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ["search.html", "_static/style.css"] + + +def setup(app): + app.add_css_file("style.css") diff --git a/mrinversion/__init__.py b/mrinversion/__init__.py index 260c070..6285450 100644 --- a/mrinversion/__init__.py +++ b/mrinversion/__init__.py @@ -1 +1 @@ -__version__ = "0.3.1" +__version__ = "0.3.1" diff --git a/mrinversion/kernel/__init__.py b/mrinversion/kernel/__init__.py index d6c5c68..0a6d164 100644 --- a/mrinversion/kernel/__init__.py +++ b/mrinversion/kernel/__init__.py @@ -1,3 +1,3 @@ -"""Initialization of the kernel module.""" -from mrinversion.kernel.relaxation import T1 # NOQA -from mrinversion.kernel.relaxation import T2 # NOQA +"""Initialization of the kernel module.""" +from mrinversion.kernel.relaxation import T1 # NOQA +from mrinversion.kernel.relaxation import T2 # NOQA diff --git a/mrinversion/linear_model/linear_inversion.py b/mrinversion/linear_model/linear_inversion.py index 374de5c..bf8a01f 100644 --- a/mrinversion/linear_model/linear_inversion.py +++ b/mrinversion/linear_model/linear_inversion.py @@ -1,46 +1,46 @@ -import numpy as np - -# import numpy.ma as ma - -__author__ = "Deepansh J. Srivastava" -__email__ = "srivastava.89@osu.edu" - - -def find_optimum_singular_value(s): - length = s.size - s2 = s**2.0 - sj = s2 / s2.sum() - T = sj * np.log10(sj) - T[np.where(np.isnan(T))] = 0 - log_n = np.log10(length) - log_nm1 = np.log10(length - 1.0) - entropy = (-1.0 / log_n) * T.sum() - - deltaEntropy = entropy - (entropy * log_n + T) / log_nm1 - - c = deltaEntropy.mean() - d = deltaEntropy.std() - - r = np.argmin(deltaEntropy - c + d) - return r - - -def TSVD(K): - U, S, VT = np.linalg.svd(K, full_matrices=False) - r = find_optimum_singular_value(S) - return U, S, VT, r - - -# standard deviation of noise remains unchanged after unitary transformation. -def reduced_subspace_kernel_and_data(U, S, VT, signal, sigma=None): - diagS = np.diag(S) - K_tilde = np.dot(diagS, VT) - s_tilde = np.dot(U.T, signal) - - projectedSignal = np.dot(U, s_tilde) - guess_solution = np.dot(np.dot(VT.T, np.diag(1 / S)), s_tilde) - - K_tilde = np.asfortranarray(K_tilde) - s_tilde = np.asfortranarray(s_tilde) - projectedSignal = np.asfortranarray(projectedSignal) - return K_tilde, s_tilde, projectedSignal, guess_solution +import numpy as np + +# import numpy.ma as ma + +__author__ = "Deepansh J. Srivastava" +__email__ = "srivastava.89@osu.edu" + + +def find_optimum_singular_value(s): + length = s.size + s2 = s**2.0 + sj = s2 / s2.sum() + T = sj * np.log10(sj) + T[np.where(np.isnan(T))] = 0 + log_n = np.log10(length) + log_nm1 = np.log10(length - 1.0) + entropy = (-1.0 / log_n) * T.sum() + + deltaEntropy = entropy - (entropy * log_n + T) / log_nm1 + + c = deltaEntropy.mean() + d = deltaEntropy.std() + + r = np.argmin(deltaEntropy - c + d) + return r + + +def TSVD(K): + U, S, VT = np.linalg.svd(K, full_matrices=False) + r = find_optimum_singular_value(S) + return U, S, VT, r + + +# standard deviation of noise remains unchanged after unitary transformation. +def reduced_subspace_kernel_and_data(U, S, VT, signal, sigma=None): + diagS = np.diag(S) + K_tilde = np.dot(diagS, VT) + s_tilde = np.dot(U.T, signal) + + projectedSignal = np.dot(U, s_tilde) + guess_solution = np.dot(np.dot(VT.T, np.diag(1 / S)), s_tilde) + + K_tilde = np.asfortranarray(K_tilde) + s_tilde = np.asfortranarray(s_tilde) + projectedSignal = np.asfortranarray(projectedSignal) + return K_tilde, s_tilde, projectedSignal, guess_solution diff --git a/mrinversion/tests/fista_intergration_test.py b/mrinversion/tests/fista_intergration_test.py index 03dbe5d..1290756 100644 --- a/mrinversion/tests/fista_intergration_test.py +++ b/mrinversion/tests/fista_intergration_test.py @@ -1,55 +1,55 @@ -"""Integration tests for FISTA implementation.""" -import csdmpy as cp -import numpy as np - -from mrinversion.kernel import relaxation -from mrinversion.linear_model import LassoFistaCV -from mrinversion.linear_model import TSVDCompression - - -def test_fista(): - domain = "https://ssnmr.org/resources/mrinversion" - filename = f"{domain}/test1_signal.csdf" - signal = cp.load(filename) - sigma = 0.0008 - - datafile = f"{domain}/test1_t2.csdf" - true_dist = cp.load(datafile) - - kernel_dimension = signal.dimensions[0] - relaxT2 = relaxation.T2( - kernel_dimension=kernel_dimension, - inverse_dimension=dict( - count=64, - minimum="1e-2 s", - maximum="1e3 s", - scale="log", - label="log (T2 / s)", - ), - ) - inverse_dimension = relaxT2.inverse_dimension - K = relaxT2.kernel(supersampling=1) - - new_system = TSVDCompression(K, signal) - compressed_K = new_system.compressed_K - compressed_s = new_system.compressed_s - - assert new_system.truncation_index == 29 - - lambdas = 10 ** (-5 + 4 * (np.arange(16) / 15)) - - f_lasso_cv = LassoFistaCV( - lambdas=lambdas, # A numpy array of lambda values. - folds=5, # The number of folds in n-folds cross-validation. - sigma=sigma, # noise standard deviation - inverse_dimension=inverse_dimension, # previously defined inverse dimensions. - ) - - f_lasso_cv.fit(K=compressed_K, s=compressed_s) - - sol = f_lasso_cv.f - assert np.argmax(sol.y[0].components[0]) == np.argmax(true_dist.y[0].components[0]) - - residuals = f_lasso_cv.residuals(K=K, s=signal) - std = residuals.std() - np.testing.assert_almost_equal(std.value, sigma, decimal=2) +"""Integration tests for FISTA implementation.""" +import csdmpy as cp +import numpy as np + +from mrinversion.kernel import relaxation +from mrinversion.linear_model import LassoFistaCV +from mrinversion.linear_model import TSVDCompression + + +def test_fista(): + domain = "https://ssnmr.org/resources/mrinversion" + filename = f"{domain}/test1_signal.csdf" + signal = cp.load(filename) + sigma = 0.0008 + + datafile = f"{domain}/test1_t2.csdf" + true_dist = cp.load(datafile) + + kernel_dimension = signal.dimensions[0] + relaxT2 = relaxation.T2( + kernel_dimension=kernel_dimension, + inverse_dimension=dict( + count=64, + minimum="1e-2 s", + maximum="1e3 s", + scale="log", + label="log (T2 / s)", + ), + ) + inverse_dimension = relaxT2.inverse_dimension + K = relaxT2.kernel(supersampling=1) + + new_system = TSVDCompression(K, signal) + compressed_K = new_system.compressed_K + compressed_s = new_system.compressed_s + + assert new_system.truncation_index == 29 + + lambdas = 10 ** (-5 + 4 * (np.arange(16) / 15)) + + f_lasso_cv = LassoFistaCV( + lambdas=lambdas, # A numpy array of lambda values. + folds=5, # The number of folds in n-folds cross-validation. + sigma=sigma, # noise standard deviation + inverse_dimension=inverse_dimension, # previously defined inverse dimensions. + ) + + f_lasso_cv.fit(K=compressed_K, s=compressed_s) + + sol = f_lasso_cv.f + assert np.argmax(sol.y[0].components[0]) == np.argmax(true_dist.y[0].components[0]) + + residuals = f_lasso_cv.residuals(K=K, s=signal) + std = residuals.std() + np.testing.assert_almost_equal(std.value, sigma, decimal=2) diff --git a/mrinversion/utils.py b/mrinversion/utils.py index 1a9e97c..79b8c52 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -1,623 +1,623 @@ -from copy import deepcopy -from itertools import combinations -from itertools import product - -import csdmpy as cp -import matplotlib.pyplot as plt -import numpy as np -from matplotlib import cm -from mpl_toolkits.mplot3d import Axes3D # lgtm [py/unused-import] - - -def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): - """Convert the three-dimensional p(iso, x, y) to p(iso, zeta, eta) tensor - distribution. - - Args - ---- - - csdm_object: CSDM - A CSDM object containing the 3D p(iso, x, y) distribution. - zeta: CSDM.Dimension - A CSDM dimension object describing the zeta dimension. - eta: CSDM.Dimension - A CSDM dimension object describing the eta dimension. - n: int - An integer used in linear interpolation of the data. The default is 5. - """ - [ - item.to("ppm", "nmr_frequency_ratio") - for item in csdm_object.x - if item.origin_offset != 0 - ] - data = csdm_object.y[0].components[0] - # print(f'data max: {data.max()}') - extra_dims = 1 - if len(csdm_object.x) > 2: - extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) - data.shape = (extra_dims, data.shape[-2], data.shape[-1]) - - reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) - dx = reg_x[1] - reg_x[0] - dy = reg_y[1] - reg_y[0] - sol = np.zeros((extra_dims, zeta.count, eta.count)) - - bins = [zeta.count, eta.count] - d_zeta = zeta.increment.value / 2 - d_eta = eta.increment.value / 2 - range_ = [ - [zeta.coordinates[0].value - d_zeta, zeta.coordinates[-1].value + d_zeta], - [eta.coordinates[0].value - d_eta, eta.coordinates[-1].value + d_eta], - ] - - avg_range_x = (np.arange(n) - (n - 1) / 2) * dx / n - avg_range_y = (np.arange(n) - (n - 1) / 2) * dy / n - for x_item in avg_range_x: - for y_item in avg_range_y: - x__ = np.abs(reg_x + x_item) - y__ = np.abs(reg_y + y_item) - x_, y_ = np.meshgrid(x__, y__) - x_ = x_.ravel() - y_ = y_.ravel() - - zeta_grid = np.sqrt(x_**2 + y_**2) - eta_grid = np.ones(zeta_grid.shape) - - index = np.where(x_ < y_) - eta_grid[index] = (4 / np.pi) * np.arctan(x_[index] / y_[index]) - - index = np.where(x_ > y_) - zeta_grid[index] *= -1 - eta_grid[index] = (4 / np.pi) * np.arctan(y_[index] / x_[index]) - - index = np.where(x_ == y_) - np.append(zeta, -zeta_grid[index]) - np.append(eta, np.ones(index[0].size)) - for i in range(extra_dims): - weight = deepcopy(data[i]).ravel() - # print(f'weight max: {weight.max()}') - # print(f'range: {range_}') - weight[index] /= 2 - np.append(weight, weight[index]) - sol_, _, _ = np.histogram2d( - zeta_grid, eta_grid, weights=weight, bins=bins, range=range_ - ) - sol[i] += sol_ - # print(f'sol max: {sol.max()}') - # print() - sol /= n * n - - del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y - csdm_new = cp.as_csdm(np.squeeze(sol)) - csdm_new.x[0] = eta - csdm_new.x[1] = zeta - if len(csdm_object.x) > 2: - csdm_new.x[2] = csdm_object.x[2] - return csdm_new - - -def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): - """Convert the three-dimensional p(iso, Cq, eta) to p(iso, xq, yq) tensor - distribution. - - Args - ---- - - csdm_object: CSDM - A CSDM object containing the 3D p(iso, Cq, eta) distribution. - xq: CSDM.Dimension - A CSDM dimension object describing the xq dimension. - yq: CSDM.Dimension - A CSDM dimension object describing the yq dimension. - n: int - An integer used in linear interpolation of the data. The default is 5. - """ - [ - item.to("ppm", "nmr_frequency_ratio") - for item in csdm_object.x - if item.origin_offset != 0 - ] - data = csdm_object.y[0].components[0] - print(data) - # csdm_object.plot() - extra_dims = 1 - if len(csdm_object.x) > 2: - extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) - data.shape = (extra_dims, data.shape[-2], data.shape[-1]) - - reg_cq, reg_eta = (csdm_object.x[i].coordinates.value for i in range(2)) - d_cq = reg_cq[1] - reg_cq[0] - d_eta = reg_eta[1] - reg_eta[0] - sol = np.zeros((extra_dims, xq.count, yq.count)) - - bins = [xq.count, yq.count] - dx = xq.increment.value / 2 - dy = yq.increment.value / 2 - range_ = [ - [xq.coordinates[0].value - dx, xq.coordinates[-1].value + dx], - [yq.coordinates[0].value - dy, yq.coordinates[-1].value + dy], - ] - # print(range_) - avg_range_cq = (np.arange(n) - (n - 1) / 2) * d_cq / n - avg_range_eta = (np.arange(n) - (n - 1) / 2) * d_eta / n - # print(avg_range_cq) - # print(avg_range_eta) - for i, cq_item in enumerate(avg_range_cq): - for j, eta_item in enumerate(avg_range_eta): - # print(i,j) - cq__ = reg_cq + cq_item - # print(np.abs(cq__).shape) - eta__ = np.abs(reg_eta + eta_item) - cq_, eta_ = np.meshgrid(cq__, eta__) - # print(cq_) - # print() - cq_ = cq_.ravel() - eta_ = eta_.ravel() - # print(cq_.shape) - # print(eta_.shape) - # print(eta_) - # print() - - theta = np.ones(eta_.shape) * 1e10 - # print(theta.shape) - - index = np.where(cq_ > 0) - # print(f'pos index: {index}') - theta[index] = np.pi / 2 * (1 - eta_[index] / 2.0) - - index = np.where(cq_ <= 0) - # print(f'neg index: {index}') - - theta[index] = np.pi / 4.0 * eta_[index] - # print(cq_) - # print(np.where(theta * 180/np.pi==1e10)) - # print(list(zip(cq_, theta * 180/np.pi))) - - xq_grid = np.abs(cq_) * np.sin(theta) - yq_grid = np.abs(cq_) * np.cos(theta) - - # index = np.arange(xq.count) - # print(f'index: {index}') - # print(data[0].ravel()) - # print(data[0].shape) - ################eta_grid = (2 / np.pi) * np.arctan(y_ / x_) - # print(f'extra dims: {extra_dims}') - for i in range(extra_dims): - weight = deepcopy(data[i]).ravel() - # weight[index] /= 2 - # np.append(weight, weight[index]) - # plt.plot(weight) - # plt.show() - # print(np.where(weight>0.1)) - # print(range_) - sol_, _, _ = np.histogram2d( - xq_grid, yq_grid, weights=weight, bins=bins, range=range_ - ) - sol[i] += sol_ - # print(sol) - # print(np.where(sol>0.1)) - # plt.imshow(sol[0,:,:]) - # plt.imshow(sol_) - # print(xq.coordinates.shape) - # print(yq.coordinates.shape) - # plt.plot(weight) - # plt.show() - - sol /= n * n - # plt.plot(weight) - # plt.show() - del xq_grid, yq_grid, index, cq_, eta_, avg_range_cq, avg_range_eta - csdm_new = cp.as_csdm(np.squeeze(sol)) - csdm_new.x[0] = yq - csdm_new.x[1] = xq - if len(csdm_object.x) > 2: - csdm_new.x[2] = csdm_object.x[2] - return csdm_new - - -def get_polar_grids(ax, ticks=None, offset=0): - """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. - - Args: - ax: Matplotlib Axes. - ticks: Tick coordinates where radial grids are drawn. The value can be a list - or a numpy array. The default value is None. - offset: The grid is drawn at an offset away from the origin. - """ - ylim = ax.get_ylim() - xlim = ax.get_xlim() - if ticks is None: - x = np.asarray(ax.get_xticks()) - inc = x[1] - x[0] - size = x.size - x = np.arange(size + 5) * inc - else: - x = np.asarray(ticks) - - lw = 0.3 - t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) - t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) - - ax.add_artist(t1) - ax.add_artist(t2) - for x_ in x: - if x_ - offset > 0: - ax.add_artist( - plt.Circle( - (0, 0), - x_ - offset, - fill=False, - color="k", - linestyle="--", - linewidth=lw, - alpha=0.5, - ) - ) - - angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 4.0) - angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) - for ang_ in angle1: - ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) - for ang_ in angle2: - ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) - ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) - ax.set_xlim(xlim) - ax.set_ylim(ylim) - - -def get_quadpolar_grids(ax, ticks=None, offset=0): - """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. - - Args: - ax: Matplotlib Axes. - ticks: Tick coordinates where radial grids are drawn. The value can be a list - or a numpy array. The default value is None. - offset: The grid is drawn at an offset away from the origin. - """ - ylim = ax.get_ylim() - xlim = ax.get_xlim() - if ticks is None: - x = np.asarray(ax.get_xticks()) - inc = x[1] - x[0] - size = x.size - x = np.arange(size + 5) * inc - else: - x = np.asarray(ticks) - - lw = 0.3 - # t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) - # t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) - - # ax.add_artist(t1) - # ax.add_artist(t2) - for x_ in x: - if x_ - offset > 0: - ax.add_artist( - plt.Circle( - (0, 0), - x_ - offset, - fill=False, - color="k", - linestyle="--", - linewidth=lw, - alpha=0.5, - ) - ) - - angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 2.0) - # angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) - for ang_ in angle1: - ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) - # for ang_ in angle2: - # ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) - # ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) - ax.set_xlim(xlim) - ax.set_ylim(ylim) - - -def plot_3d( - ax, - csdm_objects, - elev=28, - azim=-150, - x_lim=None, - y_lim=None, - z_lim=None, - cmap=cm.PiYG, - box=False, - clip_percent=0.0, - linewidth=1, - alpha=0.15, - **kwargs, -): - r"""Generate a 3D density plot with 2D contour and 1D projections. - - Args: - ax: Matplotlib Axes to render the plot. - csdm_objects: A 3D{1} CSDM object or a list of CSDM objects holding the data. - elev: (optional) The 3D view angle, elevation angle in the z plane. - azim: (optional) The 3D view angle, azimuth angle in the x-y plane. - x_lim: (optional) The x limit given as a list, [x_min, x_max]. - y_lim: (optional) The y limit given as a list, [y_min, y_max]. - z_lim: (optional) The z limit given as a list, [z_min, z_max]. - max_2d: (Optional) The normalization factor of the 2D contour projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the - maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. - max_1d: (Optional) The normalization factor of the 1D projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the - maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. - cmap: (Optional) The colormap or a list of colormaps used in rendering the - volumetric plot. The colormap list is applied to the ordered list of csdm_objects. - The same colormap is used for the 2D contour projections. For 1D plots, the first - color in the colormap scheme is used for the line color. - box: (Optional) If True, draw a box around the 3D data region. - clip_percent: (Optional) The amplitudes of the dataset below the given percent - is made transparent for the volumetric plot. - linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. - alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D - volume. - """ - csdm_object_list = csdm_objects - if not isinstance(csdm_objects, list): - csdm_object_list = [csdm_objects] - - if not isinstance(cmap, list): - cmap = [cmap] - - fraction = np.array([abs(item).sum() for item in csdm_object_list]) - fraction /= fraction.max() - - for index, item in enumerate(csdm_object_list): - plot_3d_( - ax, - item, - elev=elev, - azim=azim, - x_lim=x_lim, - y_lim=y_lim, - z_lim=z_lim, - fraction=fraction[index], - cmap=cmap[index], - box=box, - clip_percent=clip_percent, - linewidth=linewidth, - alpha=alpha, - **kwargs, - ) - - -def plot_3d_( - ax, - csdm_object, - elev=28, - azim=-150, - x_lim=None, - y_lim=None, - z_lim=None, - fraction=1.0, - cmap=cm.PiYG, - box=False, - clip_percent=0.0, - linewidth=1, - alpha=0.15, - **kwargs, -): - r"""Generate a 3D density plot with 2D contour and 1D projections. - - Args: - ax: Matplotlib Axes to render the plot. - csdm_object: A 3D{1} CSDM object holding the data. - elev: (optional) The 3D view angle, elevation angle in the z plane. - azim: (optional) The 3D view angle, azimuth angle in the x-y plane. - x_lim: (optional) The x limit given as a list, [x_min, x_max]. - y_lim: (optional) The y limit given as a list, [y_min, y_max]. - z_lim: (optional) The z limit given as a list, [z_min, z_max]. - max_2d: (Optional) The normalization factor of the 2D contour projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the - maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. - max_1d: (Optional) The normalization factor of the 1D projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the - maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. - cmap: (Optional) The colormap used in rendering the volumetric plot. The same - colormap is used for the 2D contour projections. For 1D plots, the first - color in the colormap scheme is used for the line color. - box: (Optional) If True, draw a box around the 3D data region. - clip_percent: (Optional) The amplitudes of the dataset below the given percent - is made transparent for the volumetric plot. - linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. - alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D - volume. - """ - lw = linewidth - if isinstance(csdm_object, cp.CSDM): - f = csdm_object.y[0].components[0].T - - a_, b_, c_ = (item for item in csdm_object.x) - - a = a_.coordinates.value - b = b_.coordinates.value - c = c_.coordinates.value - - xlabel = f"{a_.axis_label} - 0" - ylabel = f"{b_.axis_label} - 1" - zlabel = f"{c_.axis_label} - 2" - - else: - f = csdm_object - a = np.arange(f.shape[0]) - b = np.arange(f.shape[1]) - c = np.arange(f.shape[2]) - - xlabel = "x" - ylabel = "y" - zlabel = "z" - - delta_a = a[1] - a[0] - delta_b = b[1] - b[0] - delta_c = c[1] - c[0] - - clr = cmap - ck = cmap(0) - facecolors = cmap(f) - facecolors[:, :, :, -1] = f * alpha - index = np.where(f < clip_percent / 100) - facecolors[:, :, :, -1][index] = 0 - facecolors.shape = (f.size, 4) - - if x_lim is None: - x_lim = [a[0], a[-1]] - if y_lim is None: - y_lim = [b[0], b[-1]] - if z_lim is None: - z_lim = [c[0], c[-1]] - - offset_scalar = 50 - height_scalar = 10 - z_offset = (z_lim[1] - z_lim[0]) / offset_scalar - z_scale = np.abs(z_lim[1] - z_lim[0]) / height_scalar - offz_n = z_lim[0] - (delta_c / 2.0) - offz_1d = z_lim[1] + z_offset - - y_offset = (y_lim[1] - y_lim[0]) / offset_scalar - y_scale = np.abs(y_lim[1] - y_lim[0]) / height_scalar - offy = y_lim[1] + (delta_b / 2.0) - offy_n_1d = y_lim[0] - y_offset - - offx = x_lim[1] + (delta_a / 2.0) - x_scale = np.abs(x_lim[1] - x_lim[0]) / height_scalar - - if azim > -90 and azim <= 0: - offx = x_lim[0] - (delta_a / 2.0) - offy_n_1d = y_lim[0] - y_offset - - if azim > 0 and azim <= 90: - offy = y_lim[0] - (delta_b / 2.0) - offy_n_1d = y_lim[1] + y_offset - offx = x_lim[0] - (delta_a / 2.0) - - if azim > 90: - offx = x_lim[1] + (delta_a / 2.0) - offy_n_1d = y_lim[1] + y_offset - offy = y_lim[0] - (delta_b / 2.0) - - ax.set_proj_type("persp") - ax.view_init(elev=elev, azim=azim) - - # 2D x-y contour projection --------------------- - levels = (np.arange(20) + 1) / 20 - x1, y1 = np.meshgrid(a, b, indexing="ij") - dist_xy = f.mean(axis=2) - dist_xy *= fraction * z_scale / np.abs(dist_xy.max()) - ax.contour( - x1, - y1, - dist_xy, - zdir="z", - offset=offz_n, - cmap=clr, - levels=z_scale * levels, - linewidths=lw, - ) - - # 2D x-z contour projection - x1, y1 = np.meshgrid(a, c, indexing="ij") - dist_xz = f.mean(axis=1) - dist_xz *= fraction * y_scale / np.abs(dist_xz.max()) - ax.contour( - x1, - dist_xz, - y1, - zdir="y", - offset=offy, - cmap=clr, - levels=y_scale * levels, - linewidths=lw, - ) - - # 1D x-axis projection from 2D x-z projection - proj_x = dist_xz.mean(axis=1) - proj_x *= fraction * z_scale / np.abs(proj_x.max()) - ax.plot(a, np.sign(z_offset) * proj_x + offz_1d, offy, zdir="y", c=ck, linewidth=lw) - - # 1D z-axis projection from 2D x-z projection - proj_z = dist_xz.mean(axis=0) - proj_z *= fraction * y_scale / np.abs(proj_z.max()) - ax.plot( - np.sign(azim) * np.sign(y_offset) * proj_z + offy_n_1d, - c, - offx, - zdir="x", - c=ck, - linewidth=lw, - ) - ax.set_xlim(z_lim) - - # 2D y-z contour projection - x1, y1 = np.meshgrid(b, c, indexing="ij") - dist_yz = f.mean(axis=0) - dist_yz *= fraction * x_scale / np.abs(dist_yz.max()) - ax.contour( - dist_yz, - x1, - y1, - zdir="x", - offset=offx, - cmap=clr, - levels=x_scale * levels, - linewidths=lw, - ) - - # 1D y-axis projection - proj_y = dist_yz.mean(axis=1) - proj_y *= fraction * z_scale / np.abs(proj_y.max()) - ax.plot(b, np.sign(z_offset) * proj_y + offz_1d, offx, zdir="x", c=ck, linewidth=lw) - - ax.set_xlim(x_lim) - ax.set_ylim(y_lim) - ax.set_zlim(z_lim) - - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_zlabel(zlabel) - - x, y, z = np.meshgrid(a, b, c, indexing="ij") - - if "s" not in kwargs: - kwargs["s"] = 300 - ax.scatter(x.flat, y.flat, z.flat, marker="X", c=facecolors, **kwargs) - - # full box - r1 = [x_lim[0] - delta_a / 2, x_lim[-1] + delta_a / 2] - r2 = [y_lim[0] - delta_b / 2, y_lim[-1] + delta_b / 2] - r3 = [z_lim[-1] - delta_c / 2, z_lim[0] + delta_c / 2] - - l_box = lw - for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): - if np.sum(np.abs(s - e)) == r1[1] - r1[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - if np.sum(np.abs(s - e)) == r2[1] - r2[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - if np.sum(np.abs(s - e)) == r3[1] - r3[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - - # draw cube - if box: - r1 = [a[0] - delta_a / 2, a[-1] + delta_a / 2] - r2 = [b[0] - delta_b / 2, b[-1] + delta_b / 2] - r3 = [c[0] - delta_c / 2, c[-1] + delta_c / 2] - - for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): - if np.sum(np.abs(s - e)) == r1[1] - r1[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - if np.sum(np.abs(s - e)) == r2[1] - r2[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - if np.sum(np.abs(s - e)) == r3[1] - r3[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - - ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) +from copy import deepcopy +from itertools import combinations +from itertools import product + +import csdmpy as cp +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import cm +from mpl_toolkits.mplot3d import Axes3D # lgtm [py/unused-import] + + +def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): + """Convert the three-dimensional p(iso, x, y) to p(iso, zeta, eta) tensor + distribution. + + Args + ---- + + csdm_object: CSDM + A CSDM object containing the 3D p(iso, x, y) distribution. + zeta: CSDM.Dimension + A CSDM dimension object describing the zeta dimension. + eta: CSDM.Dimension + A CSDM dimension object describing the eta dimension. + n: int + An integer used in linear interpolation of the data. The default is 5. + """ + [ + item.to("ppm", "nmr_frequency_ratio") + for item in csdm_object.x + if item.origin_offset != 0 + ] + data = csdm_object.y[0].components[0] + # print(f'data max: {data.max()}') + extra_dims = 1 + if len(csdm_object.x) > 2: + extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) + data.shape = (extra_dims, data.shape[-2], data.shape[-1]) + + reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) + dx = reg_x[1] - reg_x[0] + dy = reg_y[1] - reg_y[0] + sol = np.zeros((extra_dims, zeta.count, eta.count)) + + bins = [zeta.count, eta.count] + d_zeta = zeta.increment.value / 2 + d_eta = eta.increment.value / 2 + range_ = [ + [zeta.coordinates[0].value - d_zeta, zeta.coordinates[-1].value + d_zeta], + [eta.coordinates[0].value - d_eta, eta.coordinates[-1].value + d_eta], + ] + + avg_range_x = (np.arange(n) - (n - 1) / 2) * dx / n + avg_range_y = (np.arange(n) - (n - 1) / 2) * dy / n + for x_item in avg_range_x: + for y_item in avg_range_y: + x__ = np.abs(reg_x + x_item) + y__ = np.abs(reg_y + y_item) + x_, y_ = np.meshgrid(x__, y__) + x_ = x_.ravel() + y_ = y_.ravel() + + zeta_grid = np.sqrt(x_**2 + y_**2) + eta_grid = np.ones(zeta_grid.shape) + + index = np.where(x_ < y_) + eta_grid[index] = (4 / np.pi) * np.arctan(x_[index] / y_[index]) + + index = np.where(x_ > y_) + zeta_grid[index] *= -1 + eta_grid[index] = (4 / np.pi) * np.arctan(y_[index] / x_[index]) + + index = np.where(x_ == y_) + np.append(zeta, -zeta_grid[index]) + np.append(eta, np.ones(index[0].size)) + for i in range(extra_dims): + weight = deepcopy(data[i]).ravel() + # print(f'weight max: {weight.max()}') + # print(f'range: {range_}') + weight[index] /= 2 + np.append(weight, weight[index]) + sol_, _, _ = np.histogram2d( + zeta_grid, eta_grid, weights=weight, bins=bins, range=range_ + ) + sol[i] += sol_ + # print(f'sol max: {sol.max()}') + # print() + sol /= n * n + + del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y + csdm_new = cp.as_csdm(np.squeeze(sol)) + csdm_new.x[0] = eta + csdm_new.x[1] = zeta + if len(csdm_object.x) > 2: + csdm_new.x[2] = csdm_object.x[2] + return csdm_new + + +def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): + """Convert the three-dimensional p(iso, Cq, eta) to p(iso, xq, yq) tensor + distribution. + + Args + ---- + + csdm_object: CSDM + A CSDM object containing the 3D p(iso, Cq, eta) distribution. + xq: CSDM.Dimension + A CSDM dimension object describing the xq dimension. + yq: CSDM.Dimension + A CSDM dimension object describing the yq dimension. + n: int + An integer used in linear interpolation of the data. The default is 5. + """ + [ + item.to("ppm", "nmr_frequency_ratio") + for item in csdm_object.x + if item.origin_offset != 0 + ] + data = csdm_object.y[0].components[0] + print(data) + # csdm_object.plot() + extra_dims = 1 + if len(csdm_object.x) > 2: + extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) + data.shape = (extra_dims, data.shape[-2], data.shape[-1]) + + reg_cq, reg_eta = (csdm_object.x[i].coordinates.value for i in range(2)) + d_cq = reg_cq[1] - reg_cq[0] + d_eta = reg_eta[1] - reg_eta[0] + sol = np.zeros((extra_dims, xq.count, yq.count)) + + bins = [xq.count, yq.count] + dx = xq.increment.value / 2 + dy = yq.increment.value / 2 + range_ = [ + [xq.coordinates[0].value - dx, xq.coordinates[-1].value + dx], + [yq.coordinates[0].value - dy, yq.coordinates[-1].value + dy], + ] + # print(range_) + avg_range_cq = (np.arange(n) - (n - 1) / 2) * d_cq / n + avg_range_eta = (np.arange(n) - (n - 1) / 2) * d_eta / n + # print(avg_range_cq) + # print(avg_range_eta) + for i, cq_item in enumerate(avg_range_cq): + for j, eta_item in enumerate(avg_range_eta): + # print(i,j) + cq__ = reg_cq + cq_item + # print(np.abs(cq__).shape) + eta__ = np.abs(reg_eta + eta_item) + cq_, eta_ = np.meshgrid(cq__, eta__) + # print(cq_) + # print() + cq_ = cq_.ravel() + eta_ = eta_.ravel() + # print(cq_.shape) + # print(eta_.shape) + # print(eta_) + # print() + + theta = np.ones(eta_.shape) * 1e10 + # print(theta.shape) + + index = np.where(cq_ > 0) + # print(f'pos index: {index}') + theta[index] = np.pi / 2 * (1 - eta_[index] / 2.0) + + index = np.where(cq_ <= 0) + # print(f'neg index: {index}') + + theta[index] = np.pi / 4.0 * eta_[index] + # print(cq_) + # print(np.where(theta * 180/np.pi==1e10)) + # print(list(zip(cq_, theta * 180/np.pi))) + + xq_grid = np.abs(cq_) * np.sin(theta) + yq_grid = np.abs(cq_) * np.cos(theta) + + # index = np.arange(xq.count) + # print(f'index: {index}') + # print(data[0].ravel()) + # print(data[0].shape) + ################eta_grid = (2 / np.pi) * np.arctan(y_ / x_) + # print(f'extra dims: {extra_dims}') + for i in range(extra_dims): + weight = deepcopy(data[i]).ravel() + # weight[index] /= 2 + # np.append(weight, weight[index]) + # plt.plot(weight) + # plt.show() + # print(np.where(weight>0.1)) + # print(range_) + sol_, _, _ = np.histogram2d( + xq_grid, yq_grid, weights=weight, bins=bins, range=range_ + ) + sol[i] += sol_ + # print(sol) + # print(np.where(sol>0.1)) + # plt.imshow(sol[0,:,:]) + # plt.imshow(sol_) + # print(xq.coordinates.shape) + # print(yq.coordinates.shape) + # plt.plot(weight) + # plt.show() + + sol /= n * n + # plt.plot(weight) + # plt.show() + del xq_grid, yq_grid, index, cq_, eta_, avg_range_cq, avg_range_eta + csdm_new = cp.as_csdm(np.squeeze(sol)) + csdm_new.x[0] = yq + csdm_new.x[1] = xq + if len(csdm_object.x) > 2: + csdm_new.x[2] = csdm_object.x[2] + return csdm_new + + +def get_polar_grids(ax, ticks=None, offset=0): + """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. + + Args: + ax: Matplotlib Axes. + ticks: Tick coordinates where radial grids are drawn. The value can be a list + or a numpy array. The default value is None. + offset: The grid is drawn at an offset away from the origin. + """ + ylim = ax.get_ylim() + xlim = ax.get_xlim() + if ticks is None: + x = np.asarray(ax.get_xticks()) + inc = x[1] - x[0] + size = x.size + x = np.arange(size + 5) * inc + else: + x = np.asarray(ticks) + + lw = 0.3 + t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) + t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) + + ax.add_artist(t1) + ax.add_artist(t2) + for x_ in x: + if x_ - offset > 0: + ax.add_artist( + plt.Circle( + (0, 0), + x_ - offset, + fill=False, + color="k", + linestyle="--", + linewidth=lw, + alpha=0.5, + ) + ) + + angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 4.0) + angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) + for ang_ in angle1: + ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) + for ang_ in angle2: + ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) + ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + +def get_quadpolar_grids(ax, ticks=None, offset=0): + """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. + + Args: + ax: Matplotlib Axes. + ticks: Tick coordinates where radial grids are drawn. The value can be a list + or a numpy array. The default value is None. + offset: The grid is drawn at an offset away from the origin. + """ + ylim = ax.get_ylim() + xlim = ax.get_xlim() + if ticks is None: + x = np.asarray(ax.get_xticks()) + inc = x[1] - x[0] + size = x.size + x = np.arange(size + 5) * inc + else: + x = np.asarray(ticks) + + lw = 0.3 + # t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) + # t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) + + # ax.add_artist(t1) + # ax.add_artist(t2) + for x_ in x: + if x_ - offset > 0: + ax.add_artist( + plt.Circle( + (0, 0), + x_ - offset, + fill=False, + color="k", + linestyle="--", + linewidth=lw, + alpha=0.5, + ) + ) + + angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 2.0) + # angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) + for ang_ in angle1: + ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) + # for ang_ in angle2: + # ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) + # ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + +def plot_3d( + ax, + csdm_objects, + elev=28, + azim=-150, + x_lim=None, + y_lim=None, + z_lim=None, + cmap=cm.PiYG, + box=False, + clip_percent=0.0, + linewidth=1, + alpha=0.15, + **kwargs, +): + r"""Generate a 3D density plot with 2D contour and 1D projections. + + Args: + ax: Matplotlib Axes to render the plot. + csdm_objects: A 3D{1} CSDM object or a list of CSDM objects holding the data. + elev: (optional) The 3D view angle, elevation angle in the z plane. + azim: (optional) The 3D view angle, azimuth angle in the x-y plane. + x_lim: (optional) The x limit given as a list, [x_min, x_max]. + y_lim: (optional) The y limit given as a list, [y_min, y_max]. + z_lim: (optional) The z limit given as a list, [z_min, z_max]. + max_2d: (Optional) The normalization factor of the 2D contour projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the + maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. + max_1d: (Optional) The normalization factor of the 1D projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the + maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. + cmap: (Optional) The colormap or a list of colormaps used in rendering the + volumetric plot. The colormap list is applied to the ordered list of csdm_objects. + The same colormap is used for the 2D contour projections. For 1D plots, the first + color in the colormap scheme is used for the line color. + box: (Optional) If True, draw a box around the 3D data region. + clip_percent: (Optional) The amplitudes of the dataset below the given percent + is made transparent for the volumetric plot. + linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. + alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D + volume. + """ + csdm_object_list = csdm_objects + if not isinstance(csdm_objects, list): + csdm_object_list = [csdm_objects] + + if not isinstance(cmap, list): + cmap = [cmap] + + fraction = np.array([abs(item).sum() for item in csdm_object_list]) + fraction /= fraction.max() + + for index, item in enumerate(csdm_object_list): + plot_3d_( + ax, + item, + elev=elev, + azim=azim, + x_lim=x_lim, + y_lim=y_lim, + z_lim=z_lim, + fraction=fraction[index], + cmap=cmap[index], + box=box, + clip_percent=clip_percent, + linewidth=linewidth, + alpha=alpha, + **kwargs, + ) + + +def plot_3d_( + ax, + csdm_object, + elev=28, + azim=-150, + x_lim=None, + y_lim=None, + z_lim=None, + fraction=1.0, + cmap=cm.PiYG, + box=False, + clip_percent=0.0, + linewidth=1, + alpha=0.15, + **kwargs, +): + r"""Generate a 3D density plot with 2D contour and 1D projections. + + Args: + ax: Matplotlib Axes to render the plot. + csdm_object: A 3D{1} CSDM object holding the data. + elev: (optional) The 3D view angle, elevation angle in the z plane. + azim: (optional) The 3D view angle, azimuth angle in the x-y plane. + x_lim: (optional) The x limit given as a list, [x_min, x_max]. + y_lim: (optional) The y limit given as a list, [y_min, y_max]. + z_lim: (optional) The z limit given as a list, [z_min, z_max]. + max_2d: (Optional) The normalization factor of the 2D contour projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the + maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. + max_1d: (Optional) The normalization factor of the 1D projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the + maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. + cmap: (Optional) The colormap used in rendering the volumetric plot. The same + colormap is used for the 2D contour projections. For 1D plots, the first + color in the colormap scheme is used for the line color. + box: (Optional) If True, draw a box around the 3D data region. + clip_percent: (Optional) The amplitudes of the dataset below the given percent + is made transparent for the volumetric plot. + linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. + alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D + volume. + """ + lw = linewidth + if isinstance(csdm_object, cp.CSDM): + f = csdm_object.y[0].components[0].T + + a_, b_, c_ = (item for item in csdm_object.x) + + a = a_.coordinates.value + b = b_.coordinates.value + c = c_.coordinates.value + + xlabel = f"{a_.axis_label} - 0" + ylabel = f"{b_.axis_label} - 1" + zlabel = f"{c_.axis_label} - 2" + + else: + f = csdm_object + a = np.arange(f.shape[0]) + b = np.arange(f.shape[1]) + c = np.arange(f.shape[2]) + + xlabel = "x" + ylabel = "y" + zlabel = "z" + + delta_a = a[1] - a[0] + delta_b = b[1] - b[0] + delta_c = c[1] - c[0] + + clr = cmap + ck = cmap(0) + facecolors = cmap(f) + facecolors[:, :, :, -1] = f * alpha + index = np.where(f < clip_percent / 100) + facecolors[:, :, :, -1][index] = 0 + facecolors.shape = (f.size, 4) + + if x_lim is None: + x_lim = [a[0], a[-1]] + if y_lim is None: + y_lim = [b[0], b[-1]] + if z_lim is None: + z_lim = [c[0], c[-1]] + + offset_scalar = 50 + height_scalar = 10 + z_offset = (z_lim[1] - z_lim[0]) / offset_scalar + z_scale = np.abs(z_lim[1] - z_lim[0]) / height_scalar + offz_n = z_lim[0] - (delta_c / 2.0) + offz_1d = z_lim[1] + z_offset + + y_offset = (y_lim[1] - y_lim[0]) / offset_scalar + y_scale = np.abs(y_lim[1] - y_lim[0]) / height_scalar + offy = y_lim[1] + (delta_b / 2.0) + offy_n_1d = y_lim[0] - y_offset + + offx = x_lim[1] + (delta_a / 2.0) + x_scale = np.abs(x_lim[1] - x_lim[0]) / height_scalar + + if azim > -90 and azim <= 0: + offx = x_lim[0] - (delta_a / 2.0) + offy_n_1d = y_lim[0] - y_offset + + if azim > 0 and azim <= 90: + offy = y_lim[0] - (delta_b / 2.0) + offy_n_1d = y_lim[1] + y_offset + offx = x_lim[0] - (delta_a / 2.0) + + if azim > 90: + offx = x_lim[1] + (delta_a / 2.0) + offy_n_1d = y_lim[1] + y_offset + offy = y_lim[0] - (delta_b / 2.0) + + ax.set_proj_type("persp") + ax.view_init(elev=elev, azim=azim) + + # 2D x-y contour projection --------------------- + levels = (np.arange(20) + 1) / 20 + x1, y1 = np.meshgrid(a, b, indexing="ij") + dist_xy = f.mean(axis=2) + dist_xy *= fraction * z_scale / np.abs(dist_xy.max()) + ax.contour( + x1, + y1, + dist_xy, + zdir="z", + offset=offz_n, + cmap=clr, + levels=z_scale * levels, + linewidths=lw, + ) + + # 2D x-z contour projection + x1, y1 = np.meshgrid(a, c, indexing="ij") + dist_xz = f.mean(axis=1) + dist_xz *= fraction * y_scale / np.abs(dist_xz.max()) + ax.contour( + x1, + dist_xz, + y1, + zdir="y", + offset=offy, + cmap=clr, + levels=y_scale * levels, + linewidths=lw, + ) + + # 1D x-axis projection from 2D x-z projection + proj_x = dist_xz.mean(axis=1) + proj_x *= fraction * z_scale / np.abs(proj_x.max()) + ax.plot(a, np.sign(z_offset) * proj_x + offz_1d, offy, zdir="y", c=ck, linewidth=lw) + + # 1D z-axis projection from 2D x-z projection + proj_z = dist_xz.mean(axis=0) + proj_z *= fraction * y_scale / np.abs(proj_z.max()) + ax.plot( + np.sign(azim) * np.sign(y_offset) * proj_z + offy_n_1d, + c, + offx, + zdir="x", + c=ck, + linewidth=lw, + ) + ax.set_xlim(z_lim) + + # 2D y-z contour projection + x1, y1 = np.meshgrid(b, c, indexing="ij") + dist_yz = f.mean(axis=0) + dist_yz *= fraction * x_scale / np.abs(dist_yz.max()) + ax.contour( + dist_yz, + x1, + y1, + zdir="x", + offset=offx, + cmap=clr, + levels=x_scale * levels, + linewidths=lw, + ) + + # 1D y-axis projection + proj_y = dist_yz.mean(axis=1) + proj_y *= fraction * z_scale / np.abs(proj_y.max()) + ax.plot(b, np.sign(z_offset) * proj_y + offz_1d, offx, zdir="x", c=ck, linewidth=lw) + + ax.set_xlim(x_lim) + ax.set_ylim(y_lim) + ax.set_zlim(z_lim) + + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_zlabel(zlabel) + + x, y, z = np.meshgrid(a, b, c, indexing="ij") + + if "s" not in kwargs: + kwargs["s"] = 300 + ax.scatter(x.flat, y.flat, z.flat, marker="X", c=facecolors, **kwargs) + + # full box + r1 = [x_lim[0] - delta_a / 2, x_lim[-1] + delta_a / 2] + r2 = [y_lim[0] - delta_b / 2, y_lim[-1] + delta_b / 2] + r3 = [z_lim[-1] - delta_c / 2, z_lim[0] + delta_c / 2] + + l_box = lw + for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): + if np.sum(np.abs(s - e)) == r1[1] - r1[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + if np.sum(np.abs(s - e)) == r2[1] - r2[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + if np.sum(np.abs(s - e)) == r3[1] - r3[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + + # draw cube + if box: + r1 = [a[0] - delta_a / 2, a[-1] + delta_a / 2] + r2 = [b[0] - delta_b / 2, b[-1] + delta_b / 2] + r3 = [c[0] - delta_c / 2, c[-1] + delta_c / 2] + + for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): + if np.sum(np.abs(s - e)) == r1[1] - r1[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + if np.sum(np.abs(s - e)) == r2[1] - r2[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + if np.sum(np.abs(s - e)) == r3[1] - r3[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + + ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) + ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) + ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) diff --git a/setup.cfg b/setup.cfg index 0fbf5f6..a567d4a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,19 +1,19 @@ -[aliases] -test=pytest - -[tool:pytest] -addopts = - --ignore=docs/conf.py - --ignore=docs/galley_examples - --ignore=docs/_build - --ignore=examples - --ignore=setup.py - --doctest-glob='docs/*.rst' - --cov='./' - --doctest-modules - -[coverage:run] -omit = - setup.py - mrinversion/utils.py - mrinversion/linear_model/fista/__init__.py +[aliases] +test=pytest + +[tool:pytest] +addopts = + --ignore=docs/conf.py + --ignore=docs/galley_examples + --ignore=docs/_build + --ignore=examples + --ignore=setup.py + --doctest-glob='docs/*.rst' + --cov='./' + --doctest-modules + +[coverage:run] +omit = + setup.py + mrinversion/utils.py + mrinversion/linear_model/fista/__init__.py diff --git a/setup.py b/setup.py index c4432a3..004b563 100644 --- a/setup.py +++ b/setup.py @@ -1,68 +1,68 @@ -"""Setup for mrinversion package.""" -from os.path import abspath -from os.path import dirname -from os.path import join -from setuptools import find_packages, setup - - -with open("mrinversion/__init__.py", encoding="utf-8") as f: - for line in f.readlines(): - if "__version__" in line: - before_keyword, keyword, after_keyword = line.partition("=") - version = after_keyword.strip()[1:-1] - -module_dir = dirname(abspath(__file__)) - -install_requires = [ - "numpy>=2.0", - "setuptools>=27.3", - "csdmpy>=0.7", - "mrsimulator>=1.0.0", - "scikit-learn>=1.5.2", - "numba>=0.61.2", -] - -setup_requires = ["setuptools>=27.3", "numpy"] -extras = {"matplotlib": ["matplotlib>=3.0"]} - - -setup( - name="mrinversion", - version=version, - description=( - "Python based statistical learning of NMR tensor and relaxation parameters " - "distribution." - ), - long_description=open(join(module_dir, "README.md"), encoding="utf-8").read(), - long_description_content_type="text/markdown", - author="Deepansh J. Srivastava", - author_email="deepansh2012@gmail.com", - python_requires=">=3.10", - url="https://github.com/DeepanshS/mrinversion/", - packages=find_packages(), - install_requires=install_requires, - setup_requires=setup_requires, - extras_require=extras, - tests_require=["pytest", "pytest-runner"], - include_package_data=True, - zip_safe=False, - license="BSD-3-Clause", - classifiers=[ - # Trove classifiers - # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers - "Intended Audience :: Science/Research", - "Intended Audience :: Developers", - "Operating System :: OS Independent", - "Development Status :: 4 - Beta", - "License :: OSI Approved :: BSD License", - "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.10", - "Programming Language :: Python :: 3.11", - "Programming Language :: Python :: 3.12", - "Programming Language :: Python :: 3.13", - "Topic :: Scientific/Engineering", - "Topic :: Education", - "Topic :: Scientific/Engineering :: Chemistry", - "Topic :: Scientific/Engineering :: Physics", - ], -) +"""Setup for mrinversion package.""" +from os.path import abspath +from os.path import dirname +from os.path import join +from setuptools import find_packages, setup + + +with open("mrinversion/__init__.py", encoding="utf-8") as f: + for line in f.readlines(): + if "__version__" in line: + before_keyword, keyword, after_keyword = line.partition("=") + version = after_keyword.strip()[1:-1] + +module_dir = dirname(abspath(__file__)) + +install_requires = [ + "numpy>=2.0", + "setuptools>=27.3", + "csdmpy>=0.7", + "mrsimulator>=1.0.0", + "scikit-learn>=1.5.2", + "numba>=0.61.2", +] + +setup_requires = ["setuptools>=27.3", "numpy"] +extras = {"matplotlib": ["matplotlib>=3.0"]} + + +setup( + name="mrinversion", + version=version, + description=( + "Python based statistical learning of NMR tensor and relaxation parameters " + "distribution." + ), + long_description=open(join(module_dir, "README.md"), encoding="utf-8").read(), + long_description_content_type="text/markdown", + author="Deepansh J. Srivastava", + author_email="deepansh2012@gmail.com", + python_requires=">=3.10", + url="https://github.com/DeepanshS/mrinversion/", + packages=find_packages(), + install_requires=install_requires, + setup_requires=setup_requires, + extras_require=extras, + tests_require=["pytest", "pytest-runner"], + include_package_data=True, + zip_safe=False, + license="BSD-3-Clause", + classifiers=[ + # Trove classifiers + # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers + "Intended Audience :: Science/Research", + "Intended Audience :: Developers", + "Operating System :: OS Independent", + "Development Status :: 4 - Beta", + "License :: OSI Approved :: BSD License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Topic :: Scientific/Engineering", + "Topic :: Education", + "Topic :: Scientific/Engineering :: Chemistry", + "Topic :: Scientific/Engineering :: Physics", + ], +) From cc66c1fb31e678b5fbee0acece8baf7fca68ea7c Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 2 Dec 2025 06:46:20 -0500 Subject: [PATCH 20/20] lint --- mrinversion/kernel/quad_aniso.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 7f9a8eb..3f85d1f 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -171,7 +171,7 @@ def __init__( self.exp_dict = exp_dict self.anisotropic_dimension = anisotropic_dimension - def kernel(self, supersampling, eta_bound=1,nQ=3): + def kernel(self, supersampling, eta_bound=1, nQ=3): import sys sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") @@ -254,7 +254,7 @@ def __init__( self.exp_dict = exp_dict self.anisotropic_dimension = anisotropic_dimension - def kernel(self, supersampling, eta_bound=1, cq_posneg=True,n_quantum=3): + def kernel(self, supersampling, eta_bound=1, cq_posneg=True, n_quantum=3): import sys sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") @@ -265,7 +265,9 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True,n_quantum=3): isotope = self.method_args["channels"][0] - Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound,calc_pos=True) + Cq, eta, abundances = self._get_zeta_eta( + supersampling, eta_bound, calc_pos=True + ) spin_systems = [ SpinSystem( @@ -277,7 +279,7 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True,n_quantum=3): obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) spec_dim = get_spectral_dimensions(obj) - operators=smsim.build_spin_matrices(nucleus=isotope, n_quantum=n_quantum) + operators = smsim.build_spin_matrices(nucleus=isotope, n_quantum=n_quantum) amp = np.asarray( [ smsim.simulate_onesite_lineshape( @@ -288,7 +290,7 @@ def kernel(self, supersampling, eta_bound=1, cq_posneg=True,n_quantum=3): contribs="c0_c4", return_array=True, distorted=True, - operators=operators + operators=operators, ) for mysys in spin_systems ]