From e368f08f61c4ee27c4f5d27952f9a606f9115f22 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 24 Aug 2020 20:24:03 -0400 Subject: [PATCH 01/11] add relaxation functions. --- examples/sideband/plot_2D_KMg0.5O4SiO2.py | 6 +- mrinversion/kernel/relaxation.py | 33 ++++++---- mrinversion/kernel/tests/other_kernel_test.py | 20 ++---- mrinversion/linear_model/_base_l1l2.py | 64 ++++++++++++------- 4 files changed, 70 insertions(+), 53 deletions(-) diff --git a/examples/sideband/plot_2D_KMg0.5O4SiO2.py b/examples/sideband/plot_2D_KMg0.5O4SiO2.py index d5f3575..8ce9113 100644 --- a/examples/sideband/plot_2D_KMg0.5O4SiO2.py +++ b/examples/sideband/plot_2D_KMg0.5O4SiO2.py @@ -196,8 +196,8 @@ def plot2D(csdm_object, **kwargs): # increase the grid resolution for your problem if desired. # setup the pre-defined range of alpha and lambda values -lambdas = 10 ** (-5.4 - 1 * (np.arange(5) / 4)) -alphas = 10 ** (-4.5 - 1.5 * (np.arange(5) / 4)) +lambdas = 10 ** (-4.6 - 1 * (np.arange(5) / 4)) +alphas = 10 ** (-4.8 - 1.5 * (np.arange(5) / 4)) # setup the smooth lasso cross-validation class s_lasso = SmoothLassoCV( @@ -397,7 +397,7 @@ def plot2D(csdm_object, **kwargs): ax.plot(f_sol_iso, "--k", label="tensor") ax.plot(Q4_region_iso, "r", label="Q4") ax.plot(Q3_region_iso, "b", label="Q3") -ax.plot(data_iso, "-k", label="MAF") +ax.plot(data_iso, "-k", label="MAT") ax.plot(data_iso - f_sol_iso - 0.1, "gray", label="residuals") ax.set_title("Isotropic projection") ax.invert_xaxis() diff --git a/mrinversion/kernel/relaxation.py b/mrinversion/kernel/relaxation.py index b180f07..360a1cf 100644 --- a/mrinversion/kernel/relaxation.py +++ b/mrinversion/kernel/relaxation.py @@ -16,29 +16,33 @@ class T2(BaseModel): kernel_dimension: A Dimension object, or an equivalent dictionary object. This dimension must represent the pure anisotropic dimension. - inverse_kernel_dimension: A list of two Dimension objects, or equivalent + inverse_dimension: A list of two Dimension objects, or equivalent dictionary objects representing the `x`-`y` coordinate grid. """ - def __init__(self, kernel_dimension, inverse_kernel_dimension): - super().__init__(kernel_dimension, inverse_kernel_dimension, 1, 1) + def __init__(self, kernel_dimension, inverse_dimension): + # unit = kernel_dimension.coordinates.to('s') + # inverse_dimension = inverse_dimension * cp.ScalarQuantity("s") + super().__init__(kernel_dimension, inverse_dimension, 1, 1) - def kernel(self, supersampling=1): + def kernel(self): """ Return the kernel of T2 decaying functions. - Args: - supersampling: An integer. Each cell is supersampled by the factor - `supersampling`. Returns: A numpy array. """ x = self.kernel_dimension.coordinates - x_inverse = _supersampled_coordinates( - self.inverse_kernel_dimension, supersampling=supersampling - ) + x_inverse = self.inverse_kernel_dimension.coordinates + # x = self.kernel_dimension.coordinates.to("s").value + # self.inverse_kernel_dimension = ( + # self.inverse_kernel_dimension / cp.ScalarQuantity("s") + # ) + # x_inverse = _supersampled_coordinates( + # self.inverse_kernel_dimension, supersampling=supersampling + # ) amp = np.exp(np.tensordot(-(1 / x_inverse), x, 0)) - return self._averaged_kernel(amp, supersampling) + return self._averaged_kernel(amp, 1) class T1(BaseModel): @@ -52,15 +56,16 @@ class T1(BaseModel): kernel_dimension: A Dimension object, or an equivalent dictionary object. This dimension must represent the pure anisotropic dimension. - inverse_kernel_dimension: A list of two Dimension objects, or equivalent + inverse_dimension: A list of two Dimension objects, or equivalent dictionary objects representing the `x`-`y` coordinate grid. """ - def __init__(self, kernel_dimension, inverse_kernel_dimension): - super().__init__(kernel_dimension, inverse_kernel_dimension, 1, 1) + def __init__(self, kernel_dimension, inverse_dimension): + super().__init__(kernel_dimension, inverse_dimension, 1, 1) def kernel(self, supersampling=1): x = self.kernel_dimension.coordinates + x_inverse = _supersampled_coordinates( self.inverse_kernel_dimension, supersampling=supersampling ) diff --git a/mrinversion/kernel/tests/other_kernel_test.py b/mrinversion/kernel/tests/other_kernel_test.py index c1e4185..f6ee134 100644 --- a/mrinversion/kernel/tests/other_kernel_test.py +++ b/mrinversion/kernel/tests/other_kernel_test.py @@ -7,34 +7,28 @@ kernel_dimension = cp.Dimension(type="linear", count=96, increment="20 ms") -inverse_kernel_dimension = cp.Dimension( +inverse_dimension = cp.Dimension( type="monotonic", coordinates=["1ms", "10ms", "100ms", "1s", "2s"] ) def test_T2_kernel(): - T2_obj = T2( - kernel_dimension=kernel_dimension, - inverse_kernel_dimension=inverse_kernel_dimension, - ) - K = T2_obj.kernel(supersampling=1) + T2_obj = T2(kernel_dimension=kernel_dimension, inverse_dimension=inverse_dimension) + K = T2_obj.kernel() x = kernel_dimension.coordinates - x_inverse = inverse_kernel_dimension.coordinates + x_inverse = inverse_dimension.coordinates amp = np.exp(np.tensordot(-x, (1 / x_inverse), 0)) amp /= amp[:, 0].sum() assert np.allclose(K, amp) def test_T1_kernel(): - T1_obj = T1( - kernel_dimension=kernel_dimension, - inverse_kernel_dimension=inverse_kernel_dimension, - ) - K = T1_obj.kernel(supersampling=1) + T1_obj = T1(kernel_dimension=kernel_dimension, inverse_dimension=inverse_dimension) + K = T1_obj.kernel() x = kernel_dimension.coordinates - x_inverse = inverse_kernel_dimension.coordinates + x_inverse = inverse_dimension.coordinates amp = 1 - np.exp(np.tensordot(-x, (1 / x_inverse), 0)) amp /= amp[:, 0].sum() assert np.allclose(K, amp) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 9b57341..ee6aa19 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -121,22 +121,13 @@ def fit(self, K, s): copy_X=True, max_iter=self.max_iterations, tol=self.tolerance, - warm_start=False, + warm_start=True, random_state=None, selection="random", positive=self.positive, ) estimator.fit(Ks, ss) - f = estimator.coef_.copy() - if s_.shape[1] > 1: - f.shape = (s_.shape[1],) + self.f_shape - f[:, :, 0] /= 2.0 - f[:, 0, :] /= 2.0 - else: - f.shape = self.f_shape - f[:, 0] /= 2.0 - f[0, :] /= 2.0 - + f = self._get_coeff_shape(estimator.coef_.copy(), s_) f *= self.scale if isinstance(s, cp.CSDM): @@ -144,13 +135,37 @@ def fit(self, K, s): if len(s.dimensions) > 1: f.dimensions[2] = s.dimensions[1] - f.dimensions[1] = self.inverse_dimension[1] - f.dimensions[0] = self.inverse_dimension[0] + + if len(self.f_shape) == 1: + f.dimensions[0] = self.inverse_dimension[0] + if len(self.f_shape) == 2: + f.dimensions[1] = self.inverse_dimension[1] + f.dimensions[0] = self.inverse_dimension[0] self.estimator = estimator self.f = f self.n_iter = estimator.n_iter_ + def _get_coeff_shape(self, f, s_): + if len(self.f_shape) == 1: + if s_.shape[1] > 1: + f.shape = (s_.shape[1],) + self.f_shape + f[:, 0] /= 2.0 + else: + f.shape = self.f_shape + f[0] /= 2.0 + + if len(self.f_shape) == 2: + if s_.shape[1] > 1: + f.shape = (s_.shape[1],) + self.f_shape + f[:, :, 0] /= 2.0 + f[:, 0, :] /= 2.0 + else: + f.shape = self.f_shape + f[:, 0] /= 2.0 + f[0, :] /= 2.0 + return f + def predict(self, K): r""" Predict the signal using the linear model. @@ -263,7 +278,7 @@ def __init__( self.inverse_dimension = inverse_dimension self.f_shape = tuple([item.count for item in inverse_dimension])[::-1] - def fit(self, K, s): + def fit(self, K, s, scoring="neg_mean_squared_error"): 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. @@ -326,7 +341,7 @@ def fit(self, K, s): tol=self.tolerance, copy_X=True, positive=self.positive, - random_state=None, + random_state=3, warm_start=True, selection="random", ) @@ -341,7 +356,7 @@ def fit(self, K, s): for alpha_ratio_ in alpha_ratio: Ks[start_index:] *= alpha_ratio_ jobs = ( - delayed(cv)(l1_array[i], Ks, ss, cv_indexes) + delayed(cv)(l1_array[i], Ks, ss, cv_indexes, scoring) for i in range(self.cv_lambdas.size) ) self.cv_map[j] = Parallel( @@ -360,12 +375,15 @@ def fit(self, K, s): # cv_map contains negated mean square errors, therefore multiply by -1. self.cv_map *= -1 # subtract the variance. - self.cv_map -= self.sigma ** 2 + if scoring == "neg_mean_squared_error": + self.cv_map -= self.sigma ** 2 + index = np.where(self.cv_map < 0) + self.cv_map[index] = np.nan - # After subtracting the variance, any negative values in the cv grid is a - # result of fitting noise. Take the absolute value of cv to avoid such - # models. - self.cv_map = np.abs(self.cv_map) + # After subtracting the variance, any negative values in the cv grid is a + # result of fitting noise. Take the absolute value of cv to avoid such + # models. + # self.cv_map = np.abs(self.cv_map) # The argmin of the minimum value is the selected model as it has the least # prediction error. @@ -452,13 +470,13 @@ def cross_validation_curve(self): return self.cv_map -def cv(l1, X, y, cv): +def cv(l1, X, y, cv, scoring="neg_mean_squared_error"): """Return the cross-validation score as negative of mean square error.""" return cross_validate( l1, X=X, y=y, - scoring="neg_mean_squared_error", # 'neg_mean_absolute_error", + scoring=scoring, # 'neg_mean_absolute_error", cv=cv, fit_params={"check_input": False}, n_jobs=1, From 5bb71f1798d81188ef5d70d94ac7620a86a82a3c Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Fri, 5 Mar 2021 14:47:33 -0500 Subject: [PATCH 02/11] update requirements: mrsimulator, scikit-learn. update hook in .pre-commit-config.yaml --- .pre-commit-config.yaml | 25 +++++++++++++++---------- docs/requirement.rst | 4 ++-- docs/requirements.txt | 2 +- mrinversion/kernel/csa_aniso.py | 2 +- requirements-dev.txt | 4 ++-- requirements.txt | 4 ++-- setup.py | 5 +++-- 7 files changed, 26 insertions(+), 20 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a9bb8fc..b067f7a 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,28 +1,28 @@ default_language_version: - python: python3.7 + python: python default_stages: [commit, push, manual] repos: - repo: https://github.com/ambv/black - rev: stable + rev: 20.8b1 hooks: - id: black name: black entry: black - language_version: python3.7 + language_version: python require_serial: true types: [python] files: \.pyi?$ - repo: https://github.com/asottile/blacken-docs - rev: v1.0.0 + rev: v1.10.0 hooks: - id: blacken-docs - additional_dependencies: [black==19.3b0] + additional_dependencies: [black] language: python - language_version: python3.7 + language_version: python - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v1.2.3 + rev: v3.4.0 hooks: - id: check-yaml language: python @@ -32,8 +32,6 @@ repos: language: python - id: check-docstring-first language: python - - id: flake8 - language: python - id: end-of-file-fixer language: python exclude: docs/notebooks/ @@ -46,8 +44,15 @@ repos: entry: debug-statement-hook language: python types: [python] + + - repo: https://gitlab.com/pycqa/flake8 + rev: "3.8.4" + hooks: + - id: flake8 + language: python + - repo: https://github.com/asottile/reorder_python_imports - rev: v1.4.0 + rev: v2.4.0 hooks: - id: reorder-python-imports name: Reorder python imports diff --git a/docs/requirement.rst b/docs/requirement.rst index 940ee83..b5d9521 100644 --- a/docs/requirement.rst +++ b/docs/requirement.rst @@ -10,9 +10,9 @@ The ``mrinversion`` library depends on the following packages: - `numpy>=1.17 `_ - `csdmpy>=0.3.1 `_ -- `mrsimulator>=0.3 `_ (for generating +- `mrsimulator>=0.5 `_ (for generating the NMR line-shape) -- `scikit-learn>=0.22.1 `_ (for statistical leaning methods) +- `scikit-learn>=0.24.1 `_ (for statistical leaning methods) **Other packages** diff --git a/docs/requirements.txt b/docs/requirements.txt index ccb73f8..dfdd94f 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -3,4 +3,4 @@ sphinx>=2.1.2 sphinx-gallery pillow matplotlib -scikit-learn>=0.22 +scikit-learn>=0.24.1 diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index f6465cd..d44de04 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -188,7 +188,7 @@ def __init__( inverse_dimension, channel, magnetic_flux_density, - "54.735 deg", + "54.7356 deg", None, None, ) diff --git a/requirements-dev.txt b/requirements-dev.txt index a4852c5..4f398be 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,8 +2,8 @@ numpy>=1.17 matplotlib>=3.0 csdmpy>=0.3.1 joblib>=0.13.2 -mrsimulator>=0.3.0a0 -scikit-learn>=0.22 +mrsimulator>=0.5 +scikit-learn>=0.24.1 # development packages pytest-runner>=5.0 diff --git a/requirements.txt b/requirements.txt index 51e2954..2c8d068 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,5 +3,5 @@ numpy>=1.17 matplotlib>=3.0 csdmpy>=0.3.1 joblib>=0.13.2 -mrsimulator>=0.3.0a0 -scikit-learn>=0.22 +mrsimulator>=0.5 +scikit-learn>=0.24.1 diff --git a/setup.py b/setup.py index c39bffe..0930a4d 100644 --- a/setup.py +++ b/setup.py @@ -18,8 +18,8 @@ "numpy>=1.17", "setuptools>=27.3", "csdmpy>=0.3.1", - "mrsimulator>=0.3.0a0", - "scikit-learn>=0.22", + "mrsimulator>=0.5", + "scikit-learn>=0.24.1", ] setup_requires = ["setuptools>=27.3"] @@ -59,6 +59,7 @@ "Programming Language :: Python :: 3.6", "Programming Language :: Python :: 3.7", "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", "Topic :: Scientific/Engineering", ], ) From d55fd5f1cbbc6c5bcac2b9ad5fdca21f7b5a84fc Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Tue, 20 Apr 2021 14:16:22 -0400 Subject: [PATCH 03/11] -increment version to 0.1.1 -add docs for SmoothLassoLS function. --- .../workflows/continuous-integration-pip.yml | 19 ++++++---- CHANGELOG | 1 + docs/api/SmoothLassoLS.rst | 13 +++++++ docs/referenceAPI.rst | 1 + mrinversion/__init__.py | 2 +- mrinversion/linear_model/smooth_lasso.py | 38 ++++++++----------- mrinversion/utils.py | 2 +- 7 files changed, 44 insertions(+), 32 deletions(-) create mode 100644 docs/api/SmoothLassoLS.rst diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 5fbda90..2b5b017 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -14,27 +14,30 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.6, 3.7, 3.8] + python-version: [3.6, 3.7, 3.8, 3.9] steps: - - uses: actions/checkout@v2 + - name: Checkout + uses: actions/checkout@v2 - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} - - name: Install mrinversion + - name: Install dependencies run: | python -m pip install --upgrade pip pip install flake8 pytest pip install -r requirements.txt pip install -r requirements-dev.txt - python setup.py develop - name: Lint with flake8 run: | # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics --exclude="docs/*" + flake8 . --count --select=C,E,F,W,N8 --show-source --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd" # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics --exclude="docs/* mrinversion/plot.py" + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=88 --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd mrinversion/utils.py" + - name: Build and install package from source + run: pip install -e . - name: Test with pytest - run: | - pytest --cov=./ --cov-report=xml + run: pytest --cov=./ --cov-report=xml + - name: Upload coverage + uses: codecov/codecov-action@v1 diff --git a/CHANGELOG b/CHANGELOG index 47fccaf..75d289a 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -4,5 +4,6 @@ v0.1.1 What's new! ''''''''''' +- Added :func:`~mrinversion.linear_model.SmoothLassoLS` function to evaluate optimum alpha and lambda values using least-squares analysis. - 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. diff --git a/docs/api/SmoothLassoLS.rst b/docs/api/SmoothLassoLS.rst new file mode 100644 index 0000000..82664a4 --- /dev/null +++ b/docs/api/SmoothLassoLS.rst @@ -0,0 +1,13 @@ +Smooth Lasso cross-validation with Least-squares +================================================ + +.. currentmodule:: mrinversion.linear_model + +.. autoclass:: SmoothLassoLS + :show-inheritance: + + .. rubric:: Methods Documentation + + .. automethod:: fit + .. automethod:: predict + .. automethod:: residuals diff --git a/docs/referenceAPI.rst b/docs/referenceAPI.rst index b25b901..6b18184 100644 --- a/docs/referenceAPI.rst +++ b/docs/referenceAPI.rst @@ -8,5 +8,6 @@ API-Reference api/shielding_kernel api/SmoothLasso api/SmoothLassoCV + api/SmoothLassoLS api/TSVDCompression api/utils diff --git a/mrinversion/__init__.py b/mrinversion/__init__.py index 81c439d..b37a8d5 100644 --- a/mrinversion/__init__.py +++ b/mrinversion/__init__.py @@ -1,3 +1,3 @@ # -*- coding: utf-8 -*- -__version__ = "0.2.0.dev0" +__version__ = "0.1.1" diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index ee39267..79a9563 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -8,8 +8,7 @@ class SmoothLasso(GeneralL2Lasso): - r""" - The linear model trained with the combined l1 and l2 priors as the regularizer. + r"""The linear model trained with the combined l1 and l2 priors as the regularizer. The method minimizes the objective function, .. math:: @@ -96,9 +95,12 @@ def __init__( class SmoothLassoCV(GeneralL2LassoCV): - r""" - The linear model trained with the combined l1 and l2 priors as the - regularizer. The method minimizes the objective function, + r"""The linear model trained with the combined l1 and l2 priors as the regularizer + with in-built cross-validation. The cross-validation is carried out using a + stratified splitting of the signal over the given range of :math:`alpha` and + :math:`lambda` values. + + The method minimizes the objective function, .. math:: \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 @@ -132,8 +134,6 @@ class SmoothLassoCV(GeneralL2LassoCV): respective dimensions, with the constraint that :math:`\prod_{i=1}^{d}n_i = n`, where :math:`d` is the total number of dimensions. - The cross-validation is carried out using a stratified splitting of the signal. - Parameters ---------- @@ -157,8 +157,6 @@ class SmoothLassoCV(GeneralL2LassoCV): The default is True. sigma: float The standard deviation of the noise in the signal. The default is 0.0. - sigma: float - The standard deviation of the noise in the signal. The default is 0.0. randomize: bool If true, the folds are created by randomly assigning the samples to each fold. If false, a stratified sampled is used to generate folds. The default is False. @@ -171,7 +169,6 @@ class SmoothLassoCV(GeneralL2LassoCV): The number of CPUs used for computation. The default is -1, that is, all available CPUs are used. - Attributes ---------- f: ndarray or CSDM object. @@ -221,9 +218,11 @@ def __init__( class SmoothLassoLS(GeneralL2LassoLS): - r""" - The linear model trained with the combined l1 and l2 priors as the - regularizer. The method minimizes the objective function, + r"""The linear model trained with the combined l1 and l2 priors as the regularizer + with in-built cross-validation. The cross-validation is carried out using a + stratified splitting of the signal using the least-squares analysis. + + The method minimizes the objective function, .. math:: \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 @@ -257,15 +256,13 @@ class SmoothLassoLS(GeneralL2LassoLS): respective dimensions, with the constraint that :math:`\prod_{i=1}^{d}n_i = n`, where :math:`d` is the total number of dimensions. - The cross-validation is carried out using a stratified splitting of the signal. - Parameters ---------- - alphas: ndarray - A list of :math:`\alpha` hyperparameters. - lambdas: ndarray - A list of :math:`\lambda` hyperparameters. + alpha: float + A starting guess for the :math:`\alpha` hyperparameters. + lambda1: float + A starting guess for the :math:`\lambda` hyperparameters. inverse_dimension: list A list of csdmpy Dimension objects representing the inverse space. folds: int @@ -282,8 +279,6 @@ class SmoothLassoLS(GeneralL2LassoLS): The default is True. sigma: float The standard deviation of the noise in the signal. The default is 0.0. - sigma: float - The standard deviation of the noise in the signal. The default is 0.0. randomize: bool If true, the folds are created by randomly assigning the samples to each fold. If false, a stratified sampled is used to generate folds. The default is False. @@ -296,7 +291,6 @@ class SmoothLassoLS(GeneralL2LassoLS): The number of CPUs used for computation. The default is -1, that is, all available CPUs are used. - Attributes ---------- f: ndarray or CSDM object. diff --git a/mrinversion/utils.py b/mrinversion/utils.py index 0b19eaa..24be9b6 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -152,7 +152,7 @@ def plot_3d( alpha=0.15, **kwargs, ): - """Generate a 3D density plot with 2D contour and 1D projections. + r"""Generate a 3D density plot with 2D contour and 1D projections. Args: ax: Matplotlib Axes to render the plot. From 7261f3b6e0ea2f577d4de53181878dba850dd5e5 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Tue, 20 Apr 2021 14:40:52 -0400 Subject: [PATCH 04/11] example updates. --- .pre-commit-config.yaml | 1 + .../auto_examples/MAF/plot_2D_0_Rb2O2p25SiO2.ipynb | 2 +- .../notebooks/auto_examples/MAF/plot_2D_1_Na2O4p7SiO2.ipynb | 2 +- .../auto_examples/MAF/plot_2D_2_Cs2Op4p72SiO2.ipynb | 2 +- .../notebooks/auto_examples/MAF/plot_2d_3_Na2O1.5SiO2.ipynb | 2 +- docs/notebooks/auto_examples/MAF/plot_2d_4_MgO.SiO2.ipynb | 2 +- docs/notebooks/auto_examples/MAF/plot_2d_5_CaO.SiO2.ipynb | 2 +- .../auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb | 2 +- .../auto_examples/synthetic/plot_1D_1_sideband.ipynb | 2 +- .../synthetic/plot_1D_2_sideband_bimodal.ipynb | 2 +- docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb | 2 +- .../auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb | 2 +- examples/MAF/plot_2D_0_Rb2O2p25SiO2.py | 6 ++---- examples/MAF/plot_2D_1_Na2O4p7SiO2.py | 6 ++---- examples/MAF/plot_2D_2_Cs2Op4p72SiO2.py | 6 ++---- examples/MAF/plot_2d_3_Na2O1.5SiO2.py | 6 ++---- examples/MAF/plot_2d_4_MgO.SiO2.py | 6 ++---- examples/MAF/plot_2d_5_CaO.SiO2.py | 6 ++---- examples/sideband/plot_2D_KMg0.5O4SiO2.py | 6 ++---- examples/synthetic/plot_1D_1_sideband.py | 4 +--- examples/synthetic/plot_1D_2_sideband_bimodal.py | 4 +--- examples/synthetic/plot_1D_3_MAF.py | 4 +--- examples/synthetic/plot_1D_4_MAF_bimodal.py | 4 +--- mrinversion/utils.py | 2 +- 24 files changed, 31 insertions(+), 52 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b067f7a..941acbf 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -61,3 +61,4 @@ repos: language: python types: [python] minimum_pre_commit_version: "0.15.0" + exclude: (examples|setup.py) diff --git a/docs/notebooks/auto_examples/MAF/plot_2D_0_Rb2O2p25SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2D_0_Rb2O2p25SiO2.ipynb index ffd0e26..db24b5f 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2D_0_Rb2O2p25SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2D_0_Rb2O2p25SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/MAF/plot_2D_1_Na2O4p7SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2D_1_Na2O4p7SiO2.ipynb index 0c43cda..c3a247a 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2D_1_Na2O4p7SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2D_1_Na2O4p7SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLasso, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/MAF/plot_2D_2_Cs2Op4p72SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2D_2_Cs2Op4p72SiO2.ipynb index ae4e648..7bc21a1 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2D_2_Cs2Op4p72SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2D_2_Cs2Op4p72SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLasso, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/MAF/plot_2d_3_Na2O1.5SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2d_3_Na2O1.5SiO2.ipynb index 8263158..a2a35ef 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2d_3_Na2O1.5SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2d_3_Na2O1.5SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/MAF/plot_2d_4_MgO.SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2d_4_MgO.SiO2.ipynb index d1bb609..4936be1 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2d_4_MgO.SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2d_4_MgO.SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/MAF/plot_2d_5_CaO.SiO2.ipynb b/docs/notebooks/auto_examples/MAF/plot_2d_5_CaO.SiO2.ipynb index c2a3910..ea06baa 100644 --- a/docs/notebooks/auto_examples/MAF/plot_2d_5_CaO.SiO2.ipynb +++ b/docs/notebooks/auto_examples/MAF/plot_2d_5_CaO.SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb b/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb index 7059368..f73fc0c 100644 --- a/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb +++ b/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import plot_3d\nfrom mrinversion.utils import to_Haeberlen_grid" + "import csdmpy as cp\nimport csdmpy.statistics as stats\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import cm\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.kernel.utils import x_y_to_zeta_eta\nfrom mrinversion.linear_model import SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import plot_3d, to_Haeberlen_grid" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb index c756c6e..32d24b7 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb index acab5f8..563b63d 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb index e50afed..d30f9c1 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb index 08f8701..48eafbe 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso\nfrom mrinversion.linear_model import SmoothLassoCV\nfrom mrinversion.linear_model import TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { diff --git a/examples/MAF/plot_2D_0_Rb2O2p25SiO2.py b/examples/MAF/plot_2D_0_Rb2O2p25SiO2.py index 6451790..f0fe49d 100644 --- a/examples/MAF/plot_2D_0_Rb2O2p25SiO2.py +++ b/examples/MAF/plot_2D_0_Rb2O2p25SiO2.py @@ -23,10 +23,8 @@ from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.kernel.utils import x_y_to_zeta_eta -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLassoCV, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 6 diff --git a/examples/MAF/plot_2D_1_Na2O4p7SiO2.py b/examples/MAF/plot_2D_1_Na2O4p7SiO2.py index 70ca352..4298b77 100644 --- a/examples/MAF/plot_2D_1_Na2O4p7SiO2.py +++ b/examples/MAF/plot_2D_1_Na2O4p7SiO2.py @@ -23,10 +23,8 @@ from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.kernel.utils import x_y_to_zeta_eta -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLasso, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 5 diff --git a/examples/MAF/plot_2D_2_Cs2Op4p72SiO2.py b/examples/MAF/plot_2D_2_Cs2Op4p72SiO2.py index e8d61ee..897b8ef 100644 --- a/examples/MAF/plot_2D_2_Cs2Op4p72SiO2.py +++ b/examples/MAF/plot_2D_2_Cs2Op4p72SiO2.py @@ -23,10 +23,8 @@ from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.kernel.utils import x_y_to_zeta_eta -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLasso, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 5 diff --git a/examples/MAF/plot_2d_3_Na2O1.5SiO2.py b/examples/MAF/plot_2d_3_Na2O1.5SiO2.py index ab8da8b..0ae2a6c 100644 --- a/examples/MAF/plot_2d_3_Na2O1.5SiO2.py +++ b/examples/MAF/plot_2d_3_Na2O1.5SiO2.py @@ -21,10 +21,8 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLasso, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 5 diff --git a/examples/MAF/plot_2d_4_MgO.SiO2.py b/examples/MAF/plot_2d_4_MgO.SiO2.py index 7374b8d..fed0fa2 100644 --- a/examples/MAF/plot_2d_4_MgO.SiO2.py +++ b/examples/MAF/plot_2d_4_MgO.SiO2.py @@ -20,10 +20,8 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLasso, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 4 diff --git a/examples/MAF/plot_2d_5_CaO.SiO2.py b/examples/MAF/plot_2d_5_CaO.SiO2.py index 0edcb8c..02beb9f 100644 --- a/examples/MAF/plot_2d_5_CaO.SiO2.py +++ b/examples/MAF/plot_2d_5_CaO.SiO2.py @@ -20,10 +20,8 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLasso, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 4 diff --git a/examples/sideband/plot_2D_KMg0.5O4SiO2.py b/examples/sideband/plot_2D_KMg0.5O4SiO2.py index 733a833..d9e9f4c 100644 --- a/examples/sideband/plot_2D_KMg0.5O4SiO2.py +++ b/examples/sideband/plot_2D_KMg0.5O4SiO2.py @@ -20,10 +20,8 @@ from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.kernel.utils import x_y_to_zeta_eta -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression -from mrinversion.utils import plot_3d -from mrinversion.utils import to_Haeberlen_grid +from mrinversion.linear_model import SmoothLassoCV, TSVDCompression +from mrinversion.utils import plot_3d, to_Haeberlen_grid # sphinx_gallery_thumbnail_number = 6 diff --git a/examples/synthetic/plot_1D_1_sideband.py b/examples/synthetic/plot_1D_1_sideband.py index 6f27719..5d2a961 100644 --- a/examples/synthetic/plot_1D_1_sideband.py +++ b/examples/synthetic/plot_1D_1_sideband.py @@ -20,9 +20,7 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression +from mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression from mrinversion.utils import get_polar_grids # Setup for the matplotlib figures diff --git a/examples/synthetic/plot_1D_2_sideband_bimodal.py b/examples/synthetic/plot_1D_2_sideband_bimodal.py index 3aafc49..eecdb97 100644 --- a/examples/synthetic/plot_1D_2_sideband_bimodal.py +++ b/examples/synthetic/plot_1D_2_sideband_bimodal.py @@ -20,9 +20,7 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression +from mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression from mrinversion.utils import get_polar_grids # Setup for the matplotlib figures diff --git a/examples/synthetic/plot_1D_3_MAF.py b/examples/synthetic/plot_1D_3_MAF.py index f57e803..7e50579 100644 --- a/examples/synthetic/plot_1D_3_MAF.py +++ b/examples/synthetic/plot_1D_3_MAF.py @@ -20,9 +20,7 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression +from mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression from mrinversion.utils import get_polar_grids # Setup for the matplotlib figures diff --git a/examples/synthetic/plot_1D_4_MAF_bimodal.py b/examples/synthetic/plot_1D_4_MAF_bimodal.py index 94de9ed..d2597ed 100644 --- a/examples/synthetic/plot_1D_4_MAF_bimodal.py +++ b/examples/synthetic/plot_1D_4_MAF_bimodal.py @@ -20,9 +20,7 @@ from pylab import rcParams from mrinversion.kernel.nmr import ShieldingPALineshape -from mrinversion.linear_model import SmoothLasso -from mrinversion.linear_model import SmoothLassoCV -from mrinversion.linear_model import TSVDCompression +from mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression from mrinversion.utils import get_polar_grids # Setup for the matplotlib figures diff --git a/mrinversion/utils.py b/mrinversion/utils.py index 24be9b6..8eb8ed0 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -26,7 +26,7 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): n: int An interger used in linear interpolation of the data. The default is 5. """ - [item.to("ppm", "nmr_frequency_ratio") for item in csdm_object.dimensions] + _ = [item.to("ppm", "nmr_frequency_ratio") for item in csdm_object.dimensions] data = csdm_object.dependent_variables[0].components[0] iso = csdm_object.dimensions[2].coordinates.value From cc7b717406ce48f63931c43ed57979ef1154e9a5 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Tue, 20 Apr 2021 22:22:55 -0400 Subject: [PATCH 05/11] update flake8 settings. --- .flake8 | 2 +- .github/workflows/continuous-integration-pip.yml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.flake8 b/.flake8 index 1933e13..67077cc 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/utils.py filename = *.pyx, *.px*, *py max-line-length = 88 max-complexity = 10 diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 2b5b017..f9850c3 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -32,9 +32,9 @@ jobs: - name: Lint with flake8 run: | # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=C,E,F,W,N8 --show-source --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd" + flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics --exclude="docs/*" # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=88 --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd mrinversion/utils.py" + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=88 --statistics --exclude="docs/* mrinversion/utils.py" - name: Build and install package from source run: pip install -e . - name: Test with pytest From 7f4d7b3d912b8f34ad626d3201c34a29bf686921 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Wed, 21 Apr 2021 15:19:17 -0400 Subject: [PATCH 06/11] improved test coverage. --- mrinversion/linear_model/_base_l1l2.py | 703 +++++++------------ mrinversion/linear_model/smooth_lasso.py | 96 +-- mrinversion/linear_model/tsvd_compression.py | 6 +- mrinversion/tests/integration_test.py | 44 +- mrinversion/utils.py | 18 +- 5 files changed, 318 insertions(+), 549 deletions(-) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index ba28e88..cbdcb6b 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -1,10 +1,16 @@ # -*- coding: utf-8 -*- from copy import deepcopy +from typing import Any +from typing import List +from typing import Literal +from typing import Union import csdmpy as cp import numpy as np from joblib import delayed from joblib import Parallel +from pydantic import BaseModel +from pydantic import PrivateAttr from scipy.optimize import least_squares from sklearn.linear_model import Lasso from sklearn.linear_model import LassoLars @@ -20,9 +26,88 @@ __email__ = "srivastava.89@osu.edu" -class GeneralL2Lasso: - r""" - The Minimizer class solves the following equation, +class Optimizer(BaseModel): + def _get_minimizer(self, alpha): + """Return the estimator for the method""" + # The factor 0.5 for alpha in the Lasso/LassoLars problem is to compensate + # 1/(2 * n_sample) factor in OLS term. + if self.method == "multi-task": + return MultiTaskLasso( + alpha=alpha / 2.0, # self.cv_lambdas[0] / 2.0, + fit_intercept=False, + normalize=False, + # precompute=True, + max_iter=self.max_iterations, + tol=self.tolerance, + copy_X=True, + # positive=self.positive, + random_state=None, + warm_start=False, + selection="random", + ) + + if self.method == "gradient_decent": + return Lasso( + alpha=alpha / 2.0, + fit_intercept=False, + normalize=False, + precompute=True, + max_iter=self.max_iterations, + tol=self.tolerance, + copy_X=True, + positive=self.positive, + random_state=None, + warm_start=False, + selection="random", + ) + + if self.method == "lars": + return LassoLars( + alpha=alpha / 2.0, + fit_intercept=False, + verbose=True, + normalize=False, + precompute="auto", + max_iter=self.max_iterations, + copy_X=True, + fit_path=False, + positive=self.positive, + jitter=None, + random_state=0, + ) + + def _pre_fit_cv(self, K, s, alpha): + s_ = _get_proper_data(s) + self._scale = s_.max().real + s_ = s_ / self._scale + + # prod = np.asarray(self.f_shape).prod() + # if K.shape[1] != prod: + # raise ValueError( + # "The product of the shape, `f_shape`, must be equal to the length of " + # f"the axis 1 of kernel, K, {K.shape[1]} != {prod}." + # ) + half = False + if "half" in self.inverse_dimension[0].application: + half = self.inverse_dimension[0].application["half"] + + Ks, ss = _get_augmented_data( + K=K, + s=s_, + alpha=s_.size * alpha, + regularizer=self.regularizer, + f_shape=self.f_shape, + half=half, + ) + + args = (self.folds, self.regularizer, self.f_shape, self.randomize, self.times) + cv_indexes = _get_cv_indexes(K, *args) + + return Ks, ss, cv_indexes + + +class GeneralL2Lasso(Optimizer): + r"""The Minimizer class solves the following equation, .. math:: {\bf f} = \underset{{\bf f}}{\text{argmin}} \left( \frac{1}{m} \| @@ -41,6 +126,7 @@ class GeneralL2Lasso: Args: alpha: Float, the hyperparameter, :math:`\alpha`. lambda1: Float, the hyperparameter, :math:`\lambda`. + hyperparameters: Dict, a python dictionary of hyperparameters. max_iterations: Interger, the maximum number of iterations allowed when solving the problem. The default value is 10000. tolerance: Float, the tolerance at which the solution is @@ -53,37 +139,38 @@ class GeneralL2Lasso: and `sparse ridge fusion`. f_shape: The shape of the solution, :math:`{\bf f}`, given as a tuple (n1, n2, ..., nd) - Attributes: """ - def __init__( - self, - alpha=1e-3, - lambda1=1e-6, - max_iterations=10000, - tolerance=1e-5, - positive=True, - regularizer=None, - inverse_dimension=None, - method="gradient_decent", - ): - - self.hyperparameters = {"lambda": lambda1, "alpha": alpha} - self.max_iterations = max_iterations - self.tolerance = tolerance - self.positive = positive - self.regularizer = regularizer - self.inverse_dimension = inverse_dimension - self.f_shape = tuple([item.count for item in inverse_dimension])[::-1] - self.method = method - - # attributes - self.f = None - self.n_iter = None + alpha: float = 1e-3 + lambda1: float = 1e-6 + hyperparameters: dict = None + inverse_dimension: List[ + Union[cp.Dimension, cp.LinearDimension, cp.MonotonicDimension] + ] = [] + max_iterations: int = 10000 + tolerance: float = 1e-5 + positive: bool = True + regularizer: Literal["smooth lasso", "sparse ridge fusion"] = "smooth lasso" + method: Literal["multi-task", "gradient_decent", "lars"] = "gradient_decent" + f: Union[cp.CSDM, np.ndarray] = None + n_iter: int = None + _scale: float = PrivateAttr(1.0) + _estimator: Any = PrivateAttr() + + def __init__(self, **kwargs): + super().__init__(**kwargs) + if self.hyperparameters is None: + self.hyperparameters = {"lambda": self.lambda1, "alpha": self.alpha} + + class Config: + arbitrary_types_allowed = True + + @property + def f_shape(self): + return tuple([item.count for item in self.inverse_dimension])[::-1] def fit(self, K, s): - r""" - Fit the model using the coordinate descent method from scikit-learn. + r"""Fit the model using the coordinate descent method from scikit-learn. Args ---- @@ -97,14 +184,9 @@ def fit(self, K, s): """ s_ = _get_proper_data(s) - # if isinstance(s, cp.CSDM): - # self.s = s - # s_ = s.dependent_variables[0].components[0].T - # else: - # s_ = s - - # if s_.ndim == 1: - # s_ = s_[:, np.newaxis] + self._scale = s_.real.max() + s_ = s_ / self._scale + # prod = np.asarray(self.f_shape).prod() # if K.shape[1] != prod: # raise ValueError( @@ -116,26 +198,18 @@ def fit(self, K, s): if "half" in self.inverse_dimension[0].application: half = self.inverse_dimension[0].application["half"] - self.scale = s_.real.max() Ks, ss = _get_augmented_data( K=K, - s=s_ / self.scale, + s=s_, alpha=s_.size * self.hyperparameters["alpha"], regularizer=self.regularizer, f_shape=self.f_shape, half=half, ) - estimator = _get_minimizer( - self.hyperparameters["lambda"], - self.method, - self.max_iterations, - self.tolerance, - self.positive, - ) - - estimator.fit(Ks, ss) - f = estimator.coef_.copy() + _estimator = self._get_minimizer(self.hyperparameters["lambda"]) + _estimator.fit(Ks, ss) + f = _estimator.coef_.copy() if half: index = self.inverse_dimension[0].application["index"] @@ -151,33 +225,32 @@ def fit(self, K, s): f[:, 0] /= 2.0 f[0, :] /= 2.0 - f *= self.scale + f *= self._scale if isinstance(s, cp.CSDM): f = self.pack_as_csdm(f, s) # f = cp.as_csdm(f) - # if len(s.dimensions) > 1: - # f.dimensions[2] = s.dimensions[1] - # f.dimensions[1] = self.inverse_dimension[1] - # f.dimensions[0] = self.inverse_dimension[0] + # if len(s.x) > 1: + # f.x[2] = s.x[1] + # f.x[1] = self.inverse_dimension[1] + # f.x[0] = self.inverse_dimension[0] - self.estimator = estimator + self._estimator = _estimator self.f = f - self.n_iter = estimator.n_iter_ + self.n_iter = _estimator.n_iter_ def pack_as_csdm(self, f, s): f = cp.as_csdm(f) - if len(s.dimensions) > 1: - f.dimensions[2] = s.dimensions[1] - f.dimensions[1] = self.inverse_dimension[1] - f.dimensions[0] = self.inverse_dimension[0] + if len(s.x) > 1: + f.x[2] = s.x[1] + f.x[1] = self.inverse_dimension[1] + f.x[0] = self.inverse_dimension[0] return f def predict(self, K): - r""" - Predict the signal using the linear model. + r"""Predict the signal using the linear model. Args ---- @@ -191,13 +264,10 @@ def predict(self, K): ndarray A numpy array of shape (m, m_count) with the predicted values """ - predict = self.estimator.predict(K) * self.scale - - return predict + return self._estimator.predict(K) * self._scale def residuals(self, K, s): - r""" - Return the residual as the difference the data and the prediced data(fit), + r"""Return the residual as the difference the data and the prediced data(fit), following .. math:: @@ -221,11 +291,8 @@ def residuals(self, K, s): is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. csdm object """ - if isinstance(s, cp.CSDM): - s_ = s.dependent_variables[0].components[0].T - else: - s_ = s - predict = np.squeeze(self.estimator.predict(K)) * self.scale + s_ = s.y[0].components[0].T if isinstance(s, cp.CSDM) else s + predict = np.squeeze(self._estimator.predict(K)) * self._scale residue = s_ - predict if not isinstance(s, cp.CSDM): @@ -240,48 +307,20 @@ def score(self, K, s, sample_weights=None): The coefficient of determination, :math:`R^2`, of the prediction. For more information, read scikit-learn documentation. """ - return self.estimator.score(K, s / self.scale, sample_weights) - - -class GeneralL2LassoLS: - def __init__( - self, - alpha=1e-6, - lambda1=1e-6, - folds=10, - max_iterations=10000, - tolerance=1e-5, - positive=True, - sigma=0.0, - regularizer=None, - randomize=False, - times=2, - verbose=False, - inverse_dimension=None, - n_jobs=-1, - method="gradient_decent", - ): - - self.method = method - self.folds = folds - - self.n_jobs = n_jobs - self.max_iterations = max_iterations - self.tolerance = tolerance - self.positive = positive - self.sigma = sigma - self.regularizer = regularizer - self.hyperparameters = {"alpha": alpha, "lambda": lambda1} - self.f = None - self.randomize = randomize - self.times = times - self.verbose = verbose - self.inverse_dimension = inverse_dimension - self.f_shape = tuple([item.count for item in inverse_dimension])[::-1] + return self._estimator.score(K, s / self._scale, sample_weights) + + +class GeneralL2LassoLS(GeneralL2Lasso): + folds: int = 10 + sigma: float = 0.0 + randomize: bool = False + times: int = 2 + verbose: bool = False + n_jobs: int = -1 + _path: Any = PrivateAttr() def fit(self, K, s, **kwargs): - r""" - Fit the model using the coordinate descent method from scikit-learn for + 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. @@ -291,54 +330,19 @@ def fit(self, K, s, **kwargs): s: A :math:`m \times m_\text{count}` signal matrix, :math:`{\bf s}` as a csdm object or a numpy array or shape (m, m_count). """ - s_ = _get_proper_data(s) - self.scale = s_.max().real - s_ = s_ / self.scale - - cv_indexes = _get_cv_indexes( - K, - self.folds, - self.regularizer, - f_shape=self.f_shape, - random=self.randomize, - times=self.times, - ) - - half = False - if "half" in self.inverse_dimension[0].application: - half = self.inverse_dimension[0].application["half"] - - ks0 = K.shape[0] - Ks, ss = _get_augmented_data( - K=K, - s=s_, - alpha=s_.size * self.hyperparameters["alpha"], - regularizer=self.regularizer, - f_shape=self.f_shape, - half=half, - ) - # start_index = K.shape[0] - - l1 = _get_minimizer( - self.hyperparameters["lambda"], - self.method, - self.max_iterations, - self.tolerance, - self.positive, - ) + Ks, ss, cv_indexes = self._pre_fit_cv(K, s, self.hyperparameters["alpha"]) - self.path = [] + self._path = [] sigma_sq = self.sigma ** 2 + ks0 = K.shape[0] + l1 = self._get_minimizer(self.hyperparameters["lambda"]) def fnc(x0): - self.path.append(x0) - # Ks_ = deepcopy(Ks) + self._path.append(x0) Ks[ks0:] *= np.sqrt(x0[0] / self.hyperparameters["alpha"]) - self.hyperparameters["alpha"] = x0[0] + alpha = self.hyperparameters["alpha"] = x0[0] l1.alpha = x0[1] / 2.0 - mse = -cv( - l1, Ks, ss, cv_indexes, alpha=self.hyperparameters["alpha"], n_jobs=-1 - ) + mse = -cv(l1, Ks, ss, cv_indexes, alpha=alpha, n_jobs=-1) mse -= sigma_sq # mse *= 1000 return np.sqrt(2 * mse) @@ -354,21 +358,10 @@ def fnc(x0): # return avg x0 = [self.hyperparameters["alpha"], self.hyperparameters["lambda"]] - # res = minimize( - # fnc, - # x0, - # bounds=[(1e-15, 1e-4), (1e-15, 1e-4)], - # method="Powell", - # **kwargs, - # ) - res = least_squares( - fnc, - x0, - # method="Powell", - bounds=[(1e-15, 1e-15), (1e-4, 1e-4)], - verbose=2, - **kwargs, - ) + + bounds = [(1e-15, 1e-15), (1e-4, 1e-4)] + # method="Powell" + res = least_squares(fnc, x0, bounds=bounds, verbose=2, **kwargs) self.hyperparameters["alpha"] = res.x[0] self.hyperparameters["lambda"] = res.x[1] @@ -377,110 +370,64 @@ def fnc(x0): # 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 - - def predict(self, K): - r""" - Predict the signal using the linear model. - - Args: - K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). + super().fit(K, s) + + +class GeneralL2LassoCV(Optimizer): + + alphas: Union[List[float], np.ndarray] = None + lambdas: Union[List[float], np.ndarray] = None + hyperparameters: dict = None + inverse_dimension: List[ + Union[cp.Dimension, cp.LinearDimension, cp.MonotonicDimension] + ] = [] + max_iterations: int = 10000 + tolerance: float = 1e-5 + positive: bool = True + regularizer: Literal["smooth lasso", "sparse ridge fusion"] = "smooth lasso" + method: Literal["multi-task", "gradient_decent", "lars"] = "gradient_decent" + + folds: int = 10 + sigma: float = 0.0 + randomize: bool = False + times: int = 2 + verbose: bool = False + n_jobs: int = -1 + + f: Union[cp.CSDM, np.ndarray] = None + n_iter: int = None + _scale: float = PrivateAttr(1.0) + _opt: Any = PrivateAttr() + _path: Any = PrivateAttr() + _cv_alphas: Any = PrivateAttr() + _cv_lambdas: Any = PrivateAttr() + _cv_map: Any = PrivateAttr() + + class Config: + arbitrary_types_allowed = True - Return: - A numpy array of shape (m, m_count) with the predicted values. - """ - return self.opt.predict(K) - - def residuals(self, K, s): - r""" - Return the residual as the difference the data and the prediced data(fit), - following - - .. math:: - \text{residuals} = {\bf s - Kf^*} - - where :math:`{\bf f^*}` is the optimum solution. - - Args: - K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). - s: A csdm object or a :math:`m \times m_\text{count}` signal matrix, - :math:`{\bf s}`. - Return: - If `s` is a csdm object, returns a csdm object with the residuals. If `s` - is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. - """ - return self.opt.residuals(K, s) - - def score(self, K, s, sample_weights=None): - """ - Return the coefficient of determination, :math:`R^2`, of the prediction. - For more information, read scikit-learn documentation. - """ - return self.opt.score(K, s, sample_weights) - - -class GeneralL2LassoCV: - def __init__( - self, - alphas=None, - lambdas=None, - folds=10, - max_iterations=10000, - tolerance=1e-5, - positive=True, - sigma=0.0, - regularizer=None, - randomize=False, - times=2, - verbose=False, - inverse_dimension=None, - n_jobs=-1, - method="gradient_decent", - ): - - if alphas is None: - self.cv_alphas = 10 ** ((np.arange(5) / 4) * 2 - 4)[::-1] - else: - self.cv_alphas = np.asarray(alphas).ravel() + @property + def f_shape(self): + return tuple([item.count for item in self.inverse_dimension])[::-1] + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self._cv_alphas = ( + 10 ** ((np.arange(5) / 4) * 2 - 4)[::-1] + if self.alphas is None + else np.asarray(self.alphas).ravel() + ) - if lambdas is None: - self.cv_lambdas = 10 ** ((np.arange(10) / 9) * 5 - 9)[::-1] - else: - self.cv_lambdas = np.asarray(lambdas).ravel() - - self.method = method - self.folds = folds - - self.n_jobs = n_jobs - self.max_iterations = max_iterations - self.tolerance = tolerance - self.positive = positive - self.sigma = sigma - self.regularizer = regularizer - self.hyperparameters = {} - self.f = None - self.randomize = randomize - self.times = times - self.verbose = verbose - self.inverse_dimension = inverse_dimension - self.f_shape = tuple([item.count for item in inverse_dimension])[::-1] + self._cv_lambdas = ( + 10 ** ((np.arange(10) / 9) * 5 - 9)[::-1] + if self.lambdas is None + else np.asarray(self.lambdas).ravel() + ) + if self.hyperparameters is None: + self.hyperparameters = {"lambda": None, "alpha": None} def fit(self, K, s, scoring="neg_mean_squared_error"): - r""" - Fit the model using the coordinate descent method from scikit-learn for + 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. @@ -491,64 +438,20 @@ def fit(self, K, s, scoring="neg_mean_squared_error"): csdm object or a numpy array or shape (m, m_count). """ - s_ = _get_proper_data(s) - - # if isinstance(s, cp.CSDM): - # self.s = s - # s_ = s.dependent_variables[0].components[0].T - # else: - # s_ = s - - # s_ = s_[:, np.newaxis] if s_.ndim == 1 else s_ - # prod = np.asarray(self.f_shape).prod() - # if K.shape[1] != prod: - # raise ValueError( - # "The product of the shape, `f_shape`, must be equal to the length of " - # f"the axis 1 of kernel, K, {K.shape[1]} != {prod}." - # ) + Ks, ss, cv_indexes = self._pre_fit_cv(K, s, self._cv_alphas[0]) - self.scale = s_.max().real - s_ = s_ / self.scale + self._cv_map = np.zeros((self._cv_alphas.size, self._cv_lambdas.size)) - cv_indexes = _get_cv_indexes( - K, - self.folds, - self.regularizer, - f_shape=self.f_shape, - random=self.randomize, - times=self.times, - ) - self.cv_map = np.zeros((self.cv_alphas.size, self.cv_lambdas.size)) - - alpha_ratio = np.ones(self.cv_alphas.size) - if self.cv_alphas.size != 1 and self.cv_alphas[0] != 0: - alpha_ratio[1:] = np.sqrt(self.cv_alphas[1:] / self.cv_alphas[:-1]) + alpha_ratio = np.ones(self._cv_alphas.size) + if self._cv_alphas.size != 1 and self._cv_alphas[0] != 0: + alpha_ratio[1:] = np.sqrt(self._cv_alphas[1:] / self._cv_alphas[:-1]) - half = False - if "half" in self.inverse_dimension[0].application: - half = self.inverse_dimension[0].application["half"] - - Ks, ss = _get_augmented_data( - K=K, - s=s_, - alpha=s_.size * self.cv_alphas[0], - regularizer=self.regularizer, - f_shape=self.f_shape, - half=half, - ) start_index = K.shape[0] - l1 = _get_minimizer( - self.cv_lambdas[0], - self.method, - self.max_iterations, - self.tolerance, - self.positive, - ) - # l1.fit(Ks, ss) - + l1 = self._get_minimizer(self._cv_lambdas[0]) l1_array = [] - for lambda_ in self.cv_lambdas: + + for lambda_ in self._cv_lambdas: l1_array.append(deepcopy(l1)) l1_array[-1].alpha = lambda_ / 2.0 @@ -558,9 +461,9 @@ def fit(self, K, s, scoring="neg_mean_squared_error"): Ks[start_index:] *= alpha_ratio_ jobs = ( delayed(cv)(l1_, Ks, ss, cv_indexes, alpha=alpha_) - for l1_, alpha_ in zip(l1_array, self.cv_alphas) + for l1_, alpha_ in zip(l1_array, self._cv_alphas) ) - self.cv_map[j] = Parallel( + self._cv_map[j] = Parallel( n_jobs=self.n_jobs, verbose=self.verbose, **{ @@ -573,28 +476,28 @@ def fit(self, K, s, scoring="neg_mean_squared_error"): )(jobs) j += 1 - # cv_map contains negated mean square errors, therefore multiply by -1. - self.cv_map *= -1 + # _cv_map contains negated mean square errors, therefore multiply by -1. + self._cv_map *= -1 # subtract the variance. if scoring == "neg_mean_squared_error": - self.cv_map -= self.sigma ** 2 - index = np.where(self.cv_map < 0) - self.cv_map[index] = np.nan + self._cv_map -= self.sigma ** 2 + index = np.where(self._cv_map < 0) + self._cv_map[index] = np.nan # After subtracting the variance, any negative values in the cv grid is a # result of fitting noise. Take the absolute value of cv to avoid such # models. - # self.cv_map = np.abs(self.cv_map) + # self._cv_map = np.abs(self._cv_map) # The argmin of the minimum value is the selected model as it has the least # prediction error. - index = np.unravel_index(self.cv_map.argmin(), self.cv_map.shape) - self.hyperparameters["alpha"] = self.cv_alphas[index[0]] - self.hyperparameters["lambda"] = self.cv_lambdas[index[1]] + index = np.unravel_index(self._cv_map.argmin(), self._cv_map.shape) + self.hyperparameters["alpha"] = self._cv_alphas[index[0]] + self.hyperparameters["lambda"] = self._cv_lambdas[index[1]] # Calculate the solution using the complete data at the optimized lambda and # alpha values - self.opt = GeneralL2Lasso( + self._opt = GeneralL2Lasso( alpha=self.hyperparameters["alpha"], lambda1=self.hyperparameters["lambda"], max_iterations=self.max_iterations, @@ -604,52 +507,26 @@ def fit(self, K, s, scoring="neg_mean_squared_error"): inverse_dimension=self.inverse_dimension, method=self.method, ) - self.opt.fit(K, s) - self.f = self.opt.f - - # 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, - # ) - # for i in range(self.folds) - # ] - - # # print(cv_indexes) - # ks0 = K.shape[0] - # for i in range(self.folds): - # lim = ks0 - len(cv_indexes[i][1]) - # # s_csdm = cp.as_csdm(s_[cv_indexes[i][0][:lim]]) - # self.opt[i].fit(K[cv_indexes[i][0][:lim]], s_[cv_indexes[i][0][:lim]]) - - # self.f = np.asarray([self.opt[i].f for i in range(self.folds)]).mean(axis=0) - # self.f_std = np.asarray([self.opt[i].f for i in range(self.folds)]).std( - # axis=0) - - # convert cv_map to csdm - self.cv_map = cp.as_csdm(np.squeeze(self.cv_map.T)) - if len(self.cv_alphas) != 1: - d0 = cp.as_dimension(-np.log10(self.cv_alphas), label="-log(α)") - self.cv_map.dimensions[0] = d0 - - if len(self.cv_lambdas) == 1: + self._opt.fit(K, s) + self.f = self._opt.f + + # convert _cv_map to csdm + self._cv_map = cp.as_csdm(np.squeeze(self._cv_map.T)) + if len(self._cv_alphas) != 1: + d0 = cp.as_dimension(-np.log10(self._cv_alphas), label="-log(α)") + self._cv_map.x[0] = d0 + + if len(self._cv_lambdas) == 1: return - d1 = cp.as_dimension(-np.log10(self.cv_lambdas), label="-log(λ)") - if len(self.cv_alphas) != 1: - self.cv_map.dimensions[1] = d1 + d1 = cp.as_dimension(-np.log10(self._cv_lambdas), label="-log(λ)") + if len(self._cv_alphas) != 1: + self._cv_map.x[1] = d1 else: - self.cv_map.dimensions[0] = d1 + self._cv_map.x[0] = d1 def predict(self, K): - r""" - Predict the signal using the linear model. + r"""Predict the signal using the linear model. Args: K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of @@ -658,11 +535,10 @@ def predict(self, K): Return: A numpy array of shape (m, m_count) with the predicted values. """ - return self.opt.predict(K) + return self._opt.predict(K) def residuals(self, K, s): - r""" - Return the residual as the difference the data and the prediced data(fit), + r"""Return the residual as the difference the data and the prediced data(fit), following .. math:: @@ -679,14 +555,14 @@ def residuals(self, K, s): If `s` is a csdm object, returns a csdm object with the residuals. If `s` is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. """ - return self.opt.residuals(K, s) + return self._opt.residuals(K, s) def score(self, K, s, sample_weights=None): """ Return the coefficient of determination, :math:`R^2`, of the prediction. For more information, read scikit-learn documentation. """ - return self.opt.score(K, s, sample_weights) + return self._opt.score(K, s, sample_weights) @property def cross_validation_curve(self): @@ -694,57 +570,7 @@ def cross_validation_curve(self): Returns: A two-dimensional CSDM object. """ - return self.cv_map - - -def _get_minimizer(alpha, method, max_iterations, tolerance, positive): - """Return the estimator for the method""" - # The factor 0.5 for alpha in the Lasso/LassoLars problem is to compensate - # 1/(2 * n_sample) factor in OLS term. - if method == "multi-task": - return MultiTaskLasso( - alpha=alpha / 2.0, # self.cv_lambdas[0] / 2.0, - fit_intercept=False, - normalize=False, - # precompute=True, - max_iter=max_iterations, - tol=tolerance, - copy_X=True, - # positive=self.positive, - random_state=None, - warm_start=False, - selection="random", - ) - - if method == "gradient_decent": - return Lasso( - alpha=alpha / 2.0, - fit_intercept=False, - normalize=False, - precompute=True, - max_iter=max_iterations, - tol=tolerance, - copy_X=True, - positive=positive, - random_state=None, - warm_start=False, - selection="random", - ) - - if method == "lars": - return LassoLars( - alpha=alpha / 2.0, - fit_intercept=False, - verbose=True, - normalize=False, - precompute="auto", - max_iter=max_iterations, - copy_X=True, - fit_path=False, - positive=positive, - jitter=None, - random_state=0, - ) + return self._cv_map def cv(l1, X, y, cv, alpha=0, n_jobs=1): @@ -814,10 +640,8 @@ def _get_cv_indexes(K, folds, regularizer, f_shape=None, random=False, times=1): """ cv_indexes = [] ks0, ks1 = K.shape - if isinstance(f_shape, int): - f_shape = (f_shape,) + f_shape = (f_shape,) if isinstance(f_shape, int) else f_shape - tr_ = [] smooth_size = _get_smooth_size(f_shape, regularizer, ks1) tr_ = (np.arange(smooth_size) + ks0).tolist() @@ -864,8 +688,7 @@ def _get_augmented_data(K, s, alpha, regularizer, f_shape=None, half=False): ks0, ks1 = K.shape ss0, ss1 = s.shape - if isinstance(f_shape, int): - f_shape = (f_shape,) + f_shape = (f_shape,) if isinstance(f_shape, int) else f_shape smooth_size = _get_smooth_size(f_shape, regularizer, ks1) @@ -939,11 +762,7 @@ def arr_j0(s0, s1): def _get_proper_data(s): - if isinstance(s, cp.CSDM): - # self.s = s - s_ = deepcopy(s.dependent_variables[0].components[0].T) - else: - s_ = deepcopy(s) - + """Extract the numpy array from the csdm, if csdm and returns a 2D array""" + s_ = deepcopy(s.y[0].components[0].T) if isinstance(s, cp.CSDM) else deepcopy(s) s_ = s_[:, np.newaxis] if s_.ndim == 1 else s_ return s_ diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index 79a9563..ec979b6 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -72,26 +72,10 @@ class SmoothLasso(GeneralL2Lasso): The number of iterations required to reach the specified tolerance. """ - def __init__( - self, - alpha, - lambda1, - inverse_dimension, - max_iterations=10000, - tolerance=1e-5, - positive=True, - method="gradient_decent", - ): - super().__init__( - alpha=alpha, - lambda1=lambda1, - max_iterations=max_iterations, - tolerance=tolerance, - positive=positive, - regularizer="smooth lasso", - inverse_dimension=inverse_dimension, - method=method, - ) + def __init__(self, **kwargs): + if "regularizer" in kwargs: + raise KeyError("regularizer is an immutable property of SmoothLasso") + super().__init__(regularizer="smooth lasso", **kwargs) class SmoothLassoCV(GeneralL2LassoCV): @@ -183,38 +167,10 @@ class SmoothLassoCV(GeneralL2LassoCV): The cross-validation error metric determined as the mean square error. """ - def __init__( - self, - alphas, - lambdas, - inverse_dimension, - folds=10, - max_iterations=10000, - tolerance=1e-5, - positive=True, - sigma=0.0, - randomize=False, - times=2, - verbose=False, - n_jobs=-1, - method="gradient_decent", - ): - super().__init__( - alphas=alphas, - lambdas=lambdas, - inverse_dimension=inverse_dimension, - folds=folds, - max_iterations=max_iterations, - tolerance=tolerance, - positive=positive, - sigma=sigma, - regularizer="smooth lasso", - randomize=randomize, - times=times, - verbose=verbose, - n_jobs=n_jobs, - method=method, - ) + def __init__(self, **kwargs): + if "regularizer" in kwargs: + raise KeyError("regularizer is an immutable property of SmoothLasso") + super().__init__(regularizer="smooth lasso", **kwargs) class SmoothLassoLS(GeneralL2LassoLS): @@ -305,35 +261,7 @@ class SmoothLassoLS(GeneralL2LassoLS): The cross-validation error metric determined as the mean square error. """ - def __init__( - self, - alpha=1e-6, - lambda1=1e-6, - inverse_dimension=None, - folds=10, - max_iterations=10000, - tolerance=1e-5, - positive=True, - sigma=0.0, - randomize=False, - times=2, - verbose=False, - n_jobs=-1, - method="gradient_decent", - ): - super().__init__( - alpha=alpha, - lambda1=lambda1, - inverse_dimension=inverse_dimension, - folds=folds, - max_iterations=max_iterations, - tolerance=tolerance, - positive=positive, - sigma=sigma, - regularizer="smooth lasso", - randomize=randomize, - times=times, - verbose=verbose, - n_jobs=n_jobs, - method=method, - ) + def __init__(self, **kwargs): + if "regularizer" in kwargs: + raise KeyError("regularizer is an immutable property of SmoothLasso") + super().__init__(regularizer="smooth lasso", **kwargs) diff --git a/mrinversion/linear_model/tsvd_compression.py b/mrinversion/linear_model/tsvd_compression.py index 77b1737..10977c6 100644 --- a/mrinversion/linear_model/tsvd_compression.py +++ b/mrinversion/linear_model/tsvd_compression.py @@ -34,7 +34,7 @@ def __init__(self, K, s, r=None): self.truncation_index = r if isinstance(s, cp.CSDM): - signal = s.dependent_variables[0].components[0].T + signal = s.y[0].components[0].T else: signal = s ( @@ -48,7 +48,7 @@ def __init__(self, K, s, r=None): if isinstance(s, cp.CSDM): self.compressed_s = cp.as_csdm(compressed_signal.T) - if len(s.dimensions) > 1: - self.compressed_s.dimensions[1] = s.dimensions[1] + if len(s.x) > 1: + self.compressed_s.x[1] = s.x[1] else: self.compressed_s = compressed_signal diff --git a/mrinversion/tests/integration_test.py b/mrinversion/tests/integration_test.py index 133e4d8..7e53ec5 100644 --- a/mrinversion/tests/integration_test.py +++ b/mrinversion/tests/integration_test.py @@ -5,20 +5,21 @@ from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.linear_model import SmoothLasso +from mrinversion.linear_model import SmoothLassoLS from mrinversion.linear_model import TSVDCompression -def test_01(): +def setup(): # The 2D MAF dataset in csdm format filename = "https://osu.box.com/shared/static/8lnwmg0dr7y6egk40c2orpkmmugh9j7c.csdf" data_object = cp.load(filename) data_object = data_object.real - _ = [item.to("ppm", "nmr_frequency_ratio") for item in data_object.dimensions] + _ = [item.to("ppm", "nmr_frequency_ratio") for item in data_object.x] data_object = data_object.T data_object_truncated = data_object[:, 155:180] - anisotropic_dimension = data_object_truncated.dimensions[0] + anisotropic_dimension = data_object_truncated.x[0] inverse_dimensions = [ cp.LinearDimension(count=25, increment="400 Hz", label="x"), cp.LinearDimension(count=25, increment="400 Hz", label="y"), @@ -41,19 +42,20 @@ def test_01(): assert new_system.truncation_index == 87 - s_lasso = SmoothLasso( - alpha=2.07e-7, lambda1=7.85e-6, inverse_dimension=inverse_dimensions - ) - s_lasso.fit(K=compressed_K, s=compressed_s) - f_sol = s_lasso.f + return inverse_dimensions, compressed_K, compressed_s, K, data_object_truncated - residuals = s_lasso.residuals(K=K, s=data_object_truncated) + +def run_test(s_lasso, K, data_object_truncated): + residuals1 = s_lasso.residuals(K=K, s=data_object_truncated) + residuals2 = s_lasso.residuals(K=K, s=data_object_truncated.y[0].components[0].T) + assert np.allclose(residuals1.y[0].components[0].T, residuals2) # assert np.allclose(residuals.mean().value, 0.00048751) - np.testing.assert_almost_equal(residuals.std().value, 0.00336372, decimal=2) + np.testing.assert_almost_equal(residuals1.std().value, 0.00336372, decimal=2) + f_sol = s_lasso.f f_sol /= f_sol.max() - [item.to("ppm", "nmr_frequency_ratio") for item in f_sol.dimensions] + _ = [item.to("ppm", "nmr_frequency_ratio") for item in f_sol.x] Q4_region = f_sol[0:8, 0:8, 3:18] Q4_region.description = "Q4 region" @@ -97,3 +99,23 @@ def test_01(): np.asarray([6.138761032132587, 7.837190479891721, 4.210912435356488]), decimal=0, ) + + +def test_01(): + inverse_dimensions, compressed_K, compressed_s, K, data_object_truncated = setup() + + s_lasso = SmoothLasso( + alpha=2.07e-7, lambda1=7.85e-6, inverse_dimension=inverse_dimensions + ) + s_lasso.fit(K=compressed_K, s=compressed_s) + run_test(s_lasso, K, data_object_truncated) + + +def test_ls(): + inverse_dimensions, compressed_K, compressed_s, K, data_object_truncated = setup() + + s_lasso_ls = SmoothLassoLS( + alpha=2.07e-7, lambda1=7.85e-6, inverse_dimension=inverse_dimensions + ) + s_lasso_ls.fit(K=compressed_K, s=compressed_s) + run_test(s_lasso_ls, K, data_object_truncated) diff --git a/mrinversion/utils.py b/mrinversion/utils.py index 8eb8ed0..4ba2683 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -26,11 +26,11 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): n: int An interger used in linear interpolation of the data. The default is 5. """ - _ = [item.to("ppm", "nmr_frequency_ratio") for item in csdm_object.dimensions] - data = csdm_object.dependent_variables[0].components[0] - iso = csdm_object.dimensions[2].coordinates.value + _ = [item.to("ppm", "nmr_frequency_ratio") for item in csdm_object.x] + data = csdm_object.y[0].components[0] + iso = csdm_object.x[2].coordinates.value - reg_x, reg_y = [csdm_object.dimensions[i].coordinates.value for i in range(2)] + 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((data.shape[0], zeta.count, eta.count)) @@ -79,9 +79,9 @@ def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y csdm_new = cp.as_csdm(sol) - csdm_new.dimensions[0] = eta - csdm_new.dimensions[1] = zeta - csdm_new.dimensions[2] = csdm_object.dimensions[2] + csdm_new.x[0] = eta + csdm_new.x[1] = zeta + csdm_new.x[2] = csdm_object.x[2] return csdm_new @@ -189,10 +189,10 @@ def plot_3d( lw = linewidth if isinstance(csdm_object, cp.CSDM): - f = csdm_object.dependent_variables[0].components[0].T + f = csdm_object.y[0].components[0].T label = csdm_object.description - a_, b_, c_ = [item for item in csdm_object.dimensions] + a_, b_, c_ = [item for item in csdm_object.x] a = a_.coordinates.value b = b_.coordinates.value From a09c4e499dea8fa919b98b2b0a3eca172b4cbd9e Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Wed, 21 Apr 2021 15:23:01 -0400 Subject: [PATCH 07/11] added pydantic to requirements.txt --- requirements-dev.txt | 1 + requirements.txt | 1 + 2 files changed, 2 insertions(+) diff --git a/requirements-dev.txt b/requirements-dev.txt index 4f398be..6958aa1 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -4,6 +4,7 @@ csdmpy>=0.3.1 joblib>=0.13.2 mrsimulator>=0.5 scikit-learn>=0.24.1 +pydantic>=1.8.1 # development packages pytest-runner>=5.0 diff --git a/requirements.txt b/requirements.txt index 2c8d068..51646ca 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,3 +5,4 @@ csdmpy>=0.3.1 joblib>=0.13.2 mrsimulator>=0.5 scikit-learn>=0.24.1 +pydantic>=1.8.1 From 6f586ef551891e3046052a3c95c19cdc77202ef6 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Wed, 21 Apr 2021 15:27:25 -0400 Subject: [PATCH 08/11] added typing_extensions to tthe requirements.txt --- mrinversion/linear_model/_base_l1l2.py | 2 +- requirements-dev.txt | 1 + requirements.txt | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index cbdcb6b..64ed3ed 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -2,7 +2,6 @@ from copy import deepcopy from typing import Any from typing import List -from typing import Literal from typing import Union import csdmpy as cp @@ -17,6 +16,7 @@ from sklearn.linear_model import MultiTaskLasso from sklearn.model_selection import cross_validate from sklearn.model_selection import KFold +from typing_extensions import Literal from mrinversion.linear_model.tsvd_compression import TSVDCompression # noqa: F401 diff --git a/requirements-dev.txt b/requirements-dev.txt index 6958aa1..4ce44eb 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -5,6 +5,7 @@ joblib>=0.13.2 mrsimulator>=0.5 scikit-learn>=0.24.1 pydantic>=1.8.1 +typing-extensions>=3.7 # development packages pytest-runner>=5.0 diff --git a/requirements.txt b/requirements.txt index 51646ca..1b125f7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,4 @@ joblib>=0.13.2 mrsimulator>=0.5 scikit-learn>=0.24.1 pydantic>=1.8.1 +typing-extensions>=3.7 From 8bfd66f0922a5240893bd3af140adf1b8646639a Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 26 Apr 2021 12:18:38 -0400 Subject: [PATCH 09/11] code organization and cleanup. --- docs/getting_started.rst | 54 +- docs/installation.rst | 4 +- docs/introduction.rst | 26 +- .../synthetic/plot_1D_1_sideband.ipynb | 2 +- .../plot_1D_2_sideband_bimodal.ipynb | 4 +- .../synthetic/plot_1D_3_MAF.ipynb | 12 +- .../synthetic/plot_1D_4_MAF_bimodal.ipynb | 2 +- examples/synthetic/plot_1D_1_sideband.py | 5 +- .../synthetic/plot_1D_2_sideband_bimodal.py | 5 +- examples/synthetic/plot_1D_3_MAF.py | 13 +- examples/synthetic/plot_1D_4_MAF_bimodal.py | 5 +- mrinversion/__init__.py | 2 +- mrinversion/linear_model/_base_l1l2.py | 646 +++++++++--------- pyproject.toml | 6 + 14 files changed, 409 insertions(+), 377 deletions(-) create mode 100644 pyproject.toml diff --git a/docs/getting_started.rst b/docs/getting_started.rst index fd07218..5dc11a5 100644 --- a/docs/getting_started.rst +++ b/docs/getting_started.rst @@ -3,7 +3,7 @@ Getting started with ``mrinversion`` ==================================== -We have put together a set of guidelines for using the `mrinversion` package. +We have put together a set of guidelines for using the *mrinversion* package. We encourage our users to follow these guidelines for consistency. Let's examine the inversion of a purely anisotropic MAS sideband spectrum into a @@ -66,10 +66,10 @@ follows, >>> filename = "https://osu.box.com/shared/static/xnlhecn8ifzcwx09f83gsh27rhc5i5l6.csdf" >>> data_object = cp.load(filename) # load the CSDM file with the csdmpy module -Here, the variable `data_object` is a `CSDM `_ +Here, the variable *data_object* is a `CSDM `_ object. The NMR spectroscopic dimension is a frequency dimension. NMR spectroscopists, however, prefer to view the spectrum on a dimensionless scale. If the -dataset dimension within the CSDM object is in frequency, you may convert it into `ppm` +dataset dimension within the CSDM object is in frequency, you may convert it into *ppm* as follows, .. plot:: @@ -80,8 +80,8 @@ as follows, >>> # convert the dimension coordinates from `Hz` to `ppm`. >>> data_object.dimensions[0].to('ppm', 'nmr_frequency_ratio') -In the above code, we convert the dimension at index 0 from `Hz` to `ppm`. For multi-dimensional -datasets, use the appropriate indexing to convert individual dimensions to `ppm`. +In the above code, we convert the dimension at index 0 from *Hz* to *ppm*. For multi-dimensional +datasets, use the appropriate indexing to convert individual dimensions to *ppm*. For comparison, let's also include the true probability distribution from which the synthetic spinning sideband dataset is derived. @@ -152,7 +152,7 @@ csdmpy library as shown below: Here, the anisotropic dimension is sampled at 625 Hz for 32 points with an offset of -10 kHz. -Similarly, we can create the CSDM dimensions needed for the `x`-`y` inversion grid as +Similarly, we can create the CSDM dimensions needed for the *x-y* inversion grid as shown below: .. plot:: @@ -166,7 +166,7 @@ shown below: ... ] Both dimensions are sampled at every 370 Hz for 25 points. The inverse dimension at -index 0 and 1 are the `x` and `y` dimensions, respectively. +index 0 and 1 are the *x* and *y* dimensions, respectively. Generating the kernel @@ -193,28 +193,28 @@ generate the kernel as follows, In the above code, the variable ``lineshapes`` is an instance of the :class:`~mrinversion.kernel.nmr.ShieldingPALineshape` class. The three required -arguments of this class are the `anisotropic_dimension`, `inverse_dimension`, and -`channel`. We have already defined the first two arguments in the previous subsection. +arguments of this class are the *anisotropic_dimension*, *inverse_dimension*, and +*channel*. We have already defined the first two arguments in the previous subsection. The value of the channel attribute is the observed nucleus. The remaining optional arguments are the metadata that describes the environment under which the spectrum is acquired. In this example, these arguments describe a :math:`^{29}\text{Si}` pure anisotropic spinning-sideband spectrum acquired at 9.4 T magnetic flux density and spinning at the magic angle (:math:`54.735^\circ`) at 625 Hz. -The value of the `rotor_frequency` argument is the effective anisotropic modulation +The value of the *rotor_frequency* argument is the effective anisotropic modulation frequency. For measurements like PASS [#f5]_, the anisotropic modulation frequency is the physical rotor frequency. For measurements like the extended chemical shift modulation sequences (XCS) [#f6]_, or its variants, where the effective anisotropic modulation frequency is lower than the physical rotor frequency, then it should be set accordingly. -The argument `number_of_sidebands` is the maximum number of sidebands that will be +The argument *number_of_sidebands* is the maximum number of sidebands that will be computed per line-shape within the kernel. For most two-dimensional isotropic vs. pure anisotropic spinning-sideband correlation spectra, the sampling along the sideband dimension is the rotor or the effective anisotropic modulation frequency. Therefore, the -`number_of_sidebands` argument is usually the number of points along the sideband +*number_of_sidebands* argument is usually the number of points along the sideband dimension. In this example, this value is 32. -Once the `ShieldingPALineshape` instance is created, use the +Once the *ShieldingPALineshape* instance is created, use the :meth:`~mrinversion.kernel.nmr.ShieldingPALineshape.kernel` method of the instance to generate the spinning sideband kernel, as follows, @@ -229,7 +229,7 @@ instance to generate the spinning sideband kernel, as follows, Here, ``K`` is the :math:`32\times 625` kernel, where the 32 is the number of samples (sideband amplitudes), and 625 is the number of features (subspectra) on the -:math:`25 \times 25` `x`-`y` grid. The argument `supersampling` is the supersampling +:math:`25 \times 25` *x-y* grid. The argument *supersampling* is the supersampling factor. In a supersampling scheme, each grid cell is averaged over a :math:`n\times n` sub-grid, where :math:`n` is the supersampling factor. A supersampling factor of 1 is equivalent to no sub-grid averaging. @@ -242,10 +242,10 @@ Often when the kernel, K, is ill-conditioned, the solution becomes unstable in the presence of the measurement noise. An ill-conditioned kernel is the one whose singular values quickly decay to zero. In such cases, we employ the truncated singular value decomposition method to approximately represent the -kernel K onto a smaller sub-space, called the `range space`, where the +kernel K onto a smaller sub-space, called the *range space*, where the sub-space kernel is relatively well-defined. We refer to this sub-space -kernel as the `compressed kernel`. Similarly, the measurement data on the -sub-space is referred to as the `compressed signal`. The compression also +kernel as the *compressed kernel*. Similarly, the measurement data on the +sub-space is referred to as the *compressed signal*. The compression also reduces the time for further computation. To compress the kernel and the data, import the :class:`~mrinversion.linear_model.TSVDCompression` class and follow, @@ -266,7 +266,7 @@ provided as the argument, when initializing the ``TSVDCompression`` class, an op truncation index is chosen using the maximum entropy method [#f7]_, which is the default behavior. The attributes :attr:`~mrinversion.linear_model.TSVDCompression.compressed_K` and :attr:`~mrinversion.linear_model.TSVDCompression.compressed_s` holds the -compressed kernel and signal, respectively. The shape of the original signal `v.s.` the +compressed kernel and signal, respectively. The shape of the original signal *v.s.* the compressed signal is .. plot:: @@ -296,9 +296,9 @@ Import the :class:`~mrinversion.linear_model.SmoothLasso` class and follow, Here, the variable ``s_lasso`` is an instance of the :class:`~mrinversion.linear_model.SmoothLasso` class. The required arguments -of this class are `alpha` and `lambda1`, corresponding to the hyperparameters +of this class are *alpha* and *lambda1*, corresponding to the hyperparameters :math:`\alpha` and :math:`\lambda`, respectively, in the Eq. :eq:`slasso`. At the -moment, we don't know the optimum value of the `alpha` and `lambda1` parameters. +moment, we don't know the optimum value of the *alpha* and *lambda1* parameters. We start with a guess value. To solve the smooth lasso problem, use the @@ -313,8 +313,8 @@ instance as follows, >>> s_lasso.fit(K=compressed_K, s=compressed_s) The two arguments of the :meth:`~mrinversion.linear_model.SmoothLasso.fit` method are -the kernel, `K`, and the signal, `s`. In the above example, we set the value of `K` as -``compressed_K``, and correspondingly the value of `s` as ``compressed_s``. You may also +the kernel, *K*, and the signal, *s*. In the above example, we set the value of *K* as +``compressed_K``, and correspondingly the value of *s* as ``compressed_s``. You may also use the uncompressed values of the kernel and signal in this method, if desired. @@ -372,7 +372,7 @@ follows, The residuals between the 1D MAS sideband spectrum and the predicted spectrum from the guess shielding tensor parameter distribution. -The argument of the `residuals` method is the kernel and the signal data. We provide the +The argument of the *residuals* method is the kernel and the signal data. We provide the original kernel, K, and signal, s, because we desire the residuals corresponding to the original data and not the compressed data. @@ -385,11 +385,11 @@ smooth LASSO estimator used here, depends on the choice of the hyperparameters. The solution shown in the above figure is when :math:`\alpha=0.01` and :math:`\lambda=1\times 10^{-4}`. Although it's a solution, it is unlikely that this is the best solution. For this, we employ the statistical learning-based model, such as the -`n`-fold cross-validation. +*n*-fold cross-validation. The :class:`~mrinversion.linear_model.SmoothLassoCV` class is designed to solve the smooth-lasso problem for a range of :math:`\alpha` and :math:`\lambda` values and -determine the best solution using the `n`-fold cross-validation. Here, we search the +determine the best solution using the *n*-fold cross-validation. Here, we search the best model on a :math:`10 \times 10` pre-defined :math:`\alpha`-:math:`\lambda` grid, using a 10-fold cross-validation statistical learning method. The :math:`\lambda` and :math:`\alpha` values are sampled uniformly on a logarithmic scale as, @@ -423,8 +423,8 @@ Setup the smooth lasso cross-validation as follows >>> s_lasso_cv.fit(K=compressed_K, s=compressed_s) The arguments of the :class:`~mrinversion.linear_model.SmoothLassoCV` is a list -of the `alpha` and `lambda` values, along with the standard deviation of the -noise, `sigma`. The value of the argument `folds` is the number of folds used in the +of the *alpha* and *lambda* values, along with the standard deviation of the +noise, *sigma*. The value of the argument *folds* is the number of folds used in the cross-validation. As before, to solve the problem, use the :meth:`~mrinversion.linear_model.SmoothLassoCV.fit` method, whose arguments are the kernel and signal. diff --git a/docs/installation.rst b/docs/installation.rst index ffbf044..7a07558 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -27,8 +27,8 @@ terminal, $ python --version -For `Mac` users, python version 3 may be installed under the name `python3`. You may replace -`python` for `python3` in the above command and all subsequent python statements. +For *Mac* users, python version 3 may be installed under the name *python3*. You may replace +*python* for *python3* in the above command and all subsequent python statements. Installing ``mrinversion`` -------------------------- diff --git a/docs/introduction.rst b/docs/introduction.rst index 7f8ecb3..042c56b 100644 --- a/docs/introduction.rst +++ b/docs/introduction.rst @@ -1,5 +1,3 @@ -.. _introduction: - ============ Introduction ============ @@ -11,12 +9,8 @@ In ``mrinversion``, we solve for the distribution of the second-rank traceless symmetric tensor principal components, through an inversion of a pure anisotropic NMR spectrum. -.. whose frequency -.. contributions are assumed to arise predominantly from the second-rank traceless -.. symmetric tensors. - In the case of the shielding tensors, the pure anisotropic frequency spectra corresponds -the cross-sections of the 2D isotropic `vs.` anisotropic correlation spectrum, such as +the cross-sections of the 2D isotropic *vs.* anisotropic correlation spectrum, such as the 2D One Pulse (TOP) MAS, phase adjusted spinning sidebands (PASS), magic-angle turning (MAT), extended chemical shift modulation (XCS), magic-angle hopping (MAH), magic-angle flipping (MAF), and Variable Angle Correlation Spectroscopy (VACSY). A key feature of all @@ -76,7 +70,7 @@ The deciding factor whether the solution :math:`{\bf f}` exists in Eq. :eq:`eq_2 is whether or not the kernel :math:`{\bf K}` is invertible. Often, most scientific problems with practical applications suffer from singular, near-singular, or ill-conditioned kernels, where :math:`{\bf K}^{-1}` doesn't exist. -Such types of problems are termed as `ill-posed`. The inversion of a purely anisotropic +Such types of problems are termed as *ill-posed*. The inversion of a purely anisotropic NMR spectrum to the distribution of the tensorial parameters is one such ill-posed problem. @@ -95,7 +89,7 @@ methods of form \|{\bf K \cdot f} - {\bf s}\|^2_2 + g({\bf f}) \right), -where :math:`\|{\bf z}\|_2` is the `l-2` norm of :math:`{\bf z}`, :math:`g({\bf f})` +where :math:`\|{\bf z}\|_2` is the *l-2* norm of :math:`{\bf z}`, :math:`g({\bf f})` is the regularization term, and :math:`{\bf f}^\dagger` is the regularized solution. The choice of the regularization term, :math:`g({\bf f})`, is often based on prior knowledge of the system for which the linear problem is defined. For anisotropic NMR @@ -154,7 +148,7 @@ respective dimensions, with the constraint that :math:`\prod_{i=1}^{d}n_i = n`, where :math:`d` is the total number of dimensions in the solution :math:`{\bf f}`, and :math:`n` is the total number of features in the kernel, :math:`{\bf K}`. -Understanding the `x-y` plot +Understanding the *x-y* plot ---------------------------- A second-rank symmetric tensor, :math:`{\bf S}`, in a three-dimensional space, is @@ -204,7 +198,7 @@ second-rank traceless tensor parameters, :math:`\zeta`-:math:`\eta`, defined as .. math:: :label: zeta_eta_def - r_\zeta = | \zeta_ | ~~~~\text{and}~~~~ + r_\zeta = | \zeta | ~~~~\text{and}~~~~ \theta = \left\{ \begin{array}{l r} \frac{\pi}{4} \eta &: \zeta \le 0, \\ \frac{\pi}{2} \left(1 - \frac{\eta}{2} \right) &: \zeta > 0. @@ -212,14 +206,14 @@ second-rank traceless tensor parameters, :math:`\zeta`-:math:`\eta`, defined as \right. Because Cartesian grids are more manageable in computation, we re-express the above polar -piece-wise grid as the `x`-`y` Cartesian grid following, +piece-wise grid as the *x-y* Cartesian grid following, .. math:: :label: x_y_def x = r_\zeta \cos\theta ~~~~\text{and}~~~~ y = r_\zeta \sin\theta. -In the `x`-`y` grid system, the basis subspectra are relatively distinguishable. The +In the *x-y* grid system, the basis subspectra are relatively distinguishable. The ``mrinversion`` library provides a utility function to render the piece-wise polar grid for your matplotlib figures. Copy-paste the following code in your script. @@ -252,7 +246,7 @@ for your matplotlib figures. Copy-paste the following code in your script. and the diagonal :math:`x=y` is :math:`\eta=1`. If you are familiar with the matplotlib library, you may notice that most code lines are -the basic matplotlib statements, except for the line that says `get_polar_grids(ax)`. +the basic matplotlib statements, except for the line that says *get_polar_grids(ax)*. The :func:`~mrinversion.utils.get_polar_grids` is a utility function that generates the piece-wise polar grid for your figures. @@ -260,7 +254,7 @@ Here, the shielding anisotropy parameter, :math:`\zeta`, is the radial dimension and the asymmetry parameter, :math:`\eta`, is the angular dimension, defined using Eqs. :eq:`zeta_eta_def` and :eq:`x_y_def`. The region in blue and red corresponds to the positive and negative values of :math:`\zeta`, where the magnitude of the anisotropy -increases radially. The `x` and the `y`-axis are :math:`\eta=0` for the negative and positive -:math:`\zeta`, respectively. When moving towards the diagonal from `x` or `y`-axes, the +increases radially. The *x* and the *y*-axis are :math:`\eta=0` for the negative and positive +:math:`\zeta`, respectively. When moving towards the diagonal from *x* or *y*-axes, the asymmetry parameter, :math:`\eta`, uniformly increase, where the diagonal is :math:`\eta=1`. diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb index 32d24b7..ba43f76 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_1_sideband.ipynb @@ -195,7 +195,7 @@ }, "outputs": [], "source": [ - "# guess alpha and lambda values.\ns_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" + "# guess alpha and lambda values.\ns_lasso = SmoothLasso(\n hyperparameters={\"alpha\": 5e-5, \"lambda\": 5e-6},\n inverse_dimension=inverse_dimension,\n)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb index 563b63d..a9a8bd4 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb @@ -177,7 +177,7 @@ }, "outputs": [], "source": [ - "# guess alpha and lambda values.\ns_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" + "# guess alpha and lambda values.\ns_lasso = SmoothLasso(\n hyperparameters={\"alpha\": 5e-5, \"lambda\": 5e-6},\n inverse_dimension=inverse_dimension,\n)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" ] }, { @@ -231,7 +231,7 @@ }, "outputs": [], "source": [ - "lambdas = 10 ** (-5 - 1 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" + "lambdas = 10 ** (-5 - 1 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n # max_iterations=50000,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s, n_repeats=2)" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb index d30f9c1..565b438 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.grid(None)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import TSVDCompression, SmoothLasso, SmoothLassoCV\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.set_xlim(0, 100)\n ax.set_ylim(0, 100)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { @@ -195,7 +195,7 @@ }, "outputs": [], "source": [ - "# guess alpha and lambda values.\ns_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" + "# guess alpha and lambda values.\ns_lasso = SmoothLasso(\n hyperparameters={\"alpha\": 5e-5, \"lambda\": 5e-6},\n inverse_dimension=inverse_dimension,\n)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" ] }, { @@ -249,7 +249,7 @@ }, "outputs": [], "source": [ - "lambdas = 10 ** (-5.2 - 1 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" + "lambdas = 10 ** (-5.2 - 1.5 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 1.5 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n max_iterations=50000,\n)\n\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" ] }, { @@ -285,14 +285,14 @@ }, "outputs": [], "source": [ - "cv_curve = s_lasso_cv.cross_validation_curve\n\n# plot of the cross-validation curve\nplt.figure(figsize=(5, 3.5))\nax = plt.subplot(projection=\"csdm\")\nax.contour(np.log10(s_lasso_cv.cross_validation_curve), levels=25)\nax.scatter(\n -np.log10(s_lasso_cv.hyperparameters[\"alpha\"]),\n -np.log10(s_lasso_cv.hyperparameters[\"lambda\"]),\n marker=\"x\",\n color=\"k\",\n)\nplt.tight_layout(pad=0.5)\nplt.show()" + "cv_curve = s_lasso_cv.cross_validation_curve\n\n# plot of the cross-validation curve\nplt.figure(figsize=(5, 3.5))\nax = plt.subplot(projection=\"csdm\")\nax.contour(np.log10(s_lasso_cv.cross_validation_curve), levels=25)\nax.scatter(\n -np.log10(s_lasso_cv.hyperparameters[\"alpha\"]),\n -np.log10(s_lasso_cv.hyperparameters[\"lambda\"]),\n marker=\"x\",\n color=\"k\",\n label=\"Grid evaluation\",\n)\nplt.legend()\nplt.tight_layout(pad=0.5)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### The optimum solution\n\nThe :attr:`~mrinversion.linear_model.SmoothLassoCV.f` attribute of the instance holds\nthe solution.\n\n" + "### The optimum solution\n\nThe :attr:`~mrinversion.linear_model.SmoothLassoLS.f` attribute of the instance holds\nthe solution.\n\n" ] }, { @@ -303,7 +303,7 @@ }, "outputs": [], "source": [ - "f_sol = s_lasso_cv.f" + "f_sol = s_lasso.f" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb index 48eafbe..a87d9d7 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_4_MAF_bimodal.ipynb @@ -177,7 +177,7 @@ }, "outputs": [], "source": [ - "# guess alpha and lambda values.\ns_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" + "# guess alpha and lambda values.\ns_lasso = SmoothLasso(\n hyperparameters={\"alpha\": 5e-5, \"lambda\": 5e-6},\n inverse_dimension=inverse_dimension,\n)\ns_lasso.fit(K=compressed_K, s=compressed_s)\nf_sol = s_lasso.f" ] }, { diff --git a/examples/synthetic/plot_1D_1_sideband.py b/examples/synthetic/plot_1D_1_sideband.py index 5d2a961..49cc7e4 100644 --- a/examples/synthetic/plot_1D_1_sideband.py +++ b/examples/synthetic/plot_1D_1_sideband.py @@ -163,7 +163,10 @@ def plot2D(ax, csdm_object, title=""): # accordingly. # guess alpha and lambda values. -s_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension) +s_lasso = SmoothLasso( + hyperparameters={"alpha": 5e-5, "lambda": 5e-6}, + inverse_dimension=inverse_dimension, +) s_lasso.fit(K=compressed_K, s=compressed_s) f_sol = s_lasso.f diff --git a/examples/synthetic/plot_1D_2_sideband_bimodal.py b/examples/synthetic/plot_1D_2_sideband_bimodal.py index eecdb97..a14b09f 100644 --- a/examples/synthetic/plot_1D_2_sideband_bimodal.py +++ b/examples/synthetic/plot_1D_2_sideband_bimodal.py @@ -146,7 +146,10 @@ def plot2D(ax, csdm_object, title=""): # accordingly. # guess alpha and lambda values. -s_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension) +s_lasso = SmoothLasso( + hyperparameters={"alpha": 5e-5, "lambda": 5e-6}, + inverse_dimension=inverse_dimension, +) s_lasso.fit(K=compressed_K, s=compressed_s) f_sol = s_lasso.f diff --git a/examples/synthetic/plot_1D_3_MAF.py b/examples/synthetic/plot_1D_3_MAF.py index 7e50579..46f9671 100644 --- a/examples/synthetic/plot_1D_3_MAF.py +++ b/examples/synthetic/plot_1D_3_MAF.py @@ -36,7 +36,8 @@ def plot2D(ax, csdm_object, title=""): levels = (np.arange(9) + 1) / 10 ax.contourf(csdm_object, cmap="gist_ncar", levels=levels) - ax.grid(None) + ax.set_xlim(0, 100) + ax.set_ylim(0, 100) ax.set_title(title) get_polar_grids(ax) ax.set_aspect("equal") @@ -162,7 +163,10 @@ def plot2D(ax, csdm_object, title=""): # accordingly. # guess alpha and lambda values. -s_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension) +s_lasso = SmoothLasso( + hyperparameters={"alpha": 5e-5, "lambda": 5e-6}, + inverse_dimension=inverse_dimension, +) s_lasso.fit(K=compressed_K, s=compressed_s) f_sol = s_lasso.f @@ -212,8 +216,8 @@ def plot2D(ax, csdm_object, title=""): # hyperparameters. # The following code generates a range of :math:`\lambda` and :math:`\alpha` values # that are uniformly sampled on the log scale. -lambdas = 10 ** (-5.2 - 1 * (np.arange(6) / 5)) -alphas = 10 ** (-4 - 2 * (np.arange(6) / 5)) +lambdas = 10 ** (-5.2 - 1.5 * (np.arange(6) / 5)) +alphas = 10 ** (-4 - 1.5 * (np.arange(6) / 5)) # set up cross validation smooth lasso method s_lasso_cv = SmoothLassoCV( @@ -255,6 +259,7 @@ def plot2D(ax, csdm_object, title=""): marker="x", color="k", ) +plt.legend() plt.tight_layout(pad=0.5) plt.show() diff --git a/examples/synthetic/plot_1D_4_MAF_bimodal.py b/examples/synthetic/plot_1D_4_MAF_bimodal.py index d2597ed..e7cbb4d 100644 --- a/examples/synthetic/plot_1D_4_MAF_bimodal.py +++ b/examples/synthetic/plot_1D_4_MAF_bimodal.py @@ -145,7 +145,10 @@ def plot2D(ax, csdm_object, title=""): # accordingly. # guess alpha and lambda values. -s_lasso = SmoothLasso(alpha=5e-5, lambda1=5e-6, inverse_dimension=inverse_dimension) +s_lasso = SmoothLasso( + hyperparameters={"alpha": 5e-5, "lambda": 5e-6}, + inverse_dimension=inverse_dimension, +) s_lasso.fit(K=compressed_K, s=compressed_s) f_sol = s_lasso.f diff --git a/mrinversion/__init__.py b/mrinversion/__init__.py index b37a8d5..5d8bf24 100644 --- a/mrinversion/__init__.py +++ b/mrinversion/__init__.py @@ -1,3 +1,3 @@ # -*- coding: utf-8 -*- -__version__ = "0.1.1" +__version__ = "0.2.0dev1" diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 64ed3ed..909db33 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -27,8 +27,90 @@ class Optimizer(BaseModel): - def _get_minimizer(self, alpha): - """Return the estimator for the method""" + inverse_dimension: List[ + Union[cp.Dimension, cp.LinearDimension, cp.MonotonicDimension] + ] + max_iterations: int = 20000 + tolerance: float = 1e-5 + positive: bool = True + regularizer: Literal["smooth lasso", "sparse ridge fusion"] = "smooth lasso" + method: Literal["multi-task", "gradient_decent", "lars"] = "gradient_decent" + + folds: int = 10 + sigma: float = 0.0 + randomize: bool = False + times: int = 2 + + f: Union[cp.CSDM, np.ndarray] = None + n_iter: int = None + + _scale: float = PrivateAttr(1.0) + _estimator: Any = PrivateAttr() + + class Config: + arbitrary_types_allowed = True + + def predict(self, K): + r"""Predict the signal using the linear model. + + Args + ---- + + K: ndarray + A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of shape + (m, n). + + Return + ------ + ndarray + A numpy array of shape (m, m_count) with the predicted values + """ + return self._estimator.predict(K) * self._scale + + def residuals(self, K, s): + r"""Return the residual as the difference the data and the prediced data(fit), + following + + .. math:: + \text{residuals} = {\bf s - Kf^*} + + where :math:`{\bf f^*}` is the optimum solution. + + Args + ---- + K: ndarray. + A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of shape + (m, n). + s: ndarray ot CSDM object. + A csdm object or a :math:`m \times m_\text{count}` signal matrix, + :math:`{\bf s}`. + + Return + ------ + ndarray or CSDM object. + If `s` is a csdm object, returns a csdm object with the residuals. If `s` + is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. + csdm object + """ + s_ = s.y[0].components[0].T if isinstance(s, cp.CSDM) else s + predict = np.squeeze(self._estimator.predict(K)) * self._scale + residue = s_ - predict + + if not isinstance(s, cp.CSDM): + return residue + + residue = cp.as_csdm(residue.T) + residue._dimensions = s._dimensions + return residue + + def score(self, K, s, sample_weights=None): + """The coefficient of determination, :math:`R^2`, of the prediction. + For more information, read scikit-learn documentation. + """ + return self._estimator.score(K, s / self._scale, sample_weights) + + def get_optimizer(self, alpha): + """Return the scikit-learn estimator for the method.""" # The factor 0.5 for alpha in the Lasso/LassoLars problem is to compensate # 1/(2 * n_sample) factor in OLS term. if self.method == "multi-task": @@ -46,37 +128,75 @@ def _get_minimizer(self, alpha): selection="random", ) + kwargs = dict( + alpha=alpha / 2.0, + fit_intercept=False, + normalize=False, + max_iter=self.max_iterations, + copy_X=True, + positive=self.positive, + random_state=45, + ) + if self.method == "gradient_decent": return Lasso( - alpha=alpha / 2.0, - fit_intercept=False, - normalize=False, precompute=True, - max_iter=self.max_iterations, tol=self.tolerance, - copy_X=True, - positive=self.positive, - random_state=None, warm_start=False, selection="random", + **kwargs, ) if self.method == "lars": return LassoLars( - alpha=alpha / 2.0, - fit_intercept=False, - verbose=True, - normalize=False, precompute="auto", - max_iter=self.max_iterations, - copy_X=True, fit_path=False, - positive=self.positive, - jitter=None, - random_state=0, + jitter=self.sigma, + verbose=True, + **kwargs, ) - def _pre_fit_cv(self, K, s, alpha): + @property + def f_shape(self): + """Shape of solution grid""" + return tuple([item.count for item in self.inverse_dimension])[::-1] + + def pack_solution_as_csdm(self, f, s): + """Pack the solution as a CSDM object. + + Args: + np.ndarray f: Solution numpy array. + cp.csdm s: CSDM dataset + """ + f = cp.as_csdm(f) + + if len(s.x) > 1: + f.x[2] = s.x[1] + f.x[1] = self.inverse_dimension[1] + f.x[0] = self.inverse_dimension[0] + return f + + def shape_solution(self, f, s_, half=False): + if half: + index = self.inverse_dimension[0].application["index"] + f_new = np.zeros((s_.shape[1], np.prod(self.f_shape)), dtype=float) + f_new[:, index] = f + f = f_new + if s_.shape[1] > 1: + f.shape = (s_.shape[1],) + self.f_shape + f[:, :, 0] /= 2.0 + f[:, 0, :] /= 2.0 + else: + f.shape = self.f_shape + f[:, 0] /= 2.0 + f[0, :] /= 2.0 + + f *= self._scale + return f + + def _fit_setup(self, K, s, alpha, get_cv_indexes=False): + """Compress the kernel and dataset, and optionally calculate the + cross-validation indexes (test-train sets).""" s_ = _get_proper_data(s) self._scale = s_.max().real s_ = s_ / self._scale @@ -91,20 +211,33 @@ def _pre_fit_cv(self, K, s, alpha): if "half" in self.inverse_dimension[0].application: half = self.inverse_dimension[0].application["half"] - Ks, ss = _get_augmented_data( - K=K, - s=s_, - alpha=s_.size * alpha, - regularizer=self.regularizer, - f_shape=self.f_shape, - half=half, - ) + args = {"regularizer": self.regularizer, "f_shape": self.f_shape, "half": half} + Ks, ss = _get_augmented_data(K=K, s=s_, alpha=s_.size * alpha, **args) + + if not get_cv_indexes: + return Ks, ss, (half, s_) args = (self.folds, self.regularizer, self.f_shape, self.randomize, self.times) cv_indexes = _get_cv_indexes(K, *args) return Ks, ss, cv_indexes + def _get_bootstrap_cv_indexes(self, cv_indexes, n_repeats, arr_length): + if n_repeats == 1: + return [cv_indexes] + + a = np.arange(arr_length) + np.random.shuffle(a) + list_ = a[:n_repeats] + + return [ + [ + [np.delete(t_index, np.where(t_index == list_[i])) for t_index in items] + for items in cv_indexes + ] + for i in range(n_repeats) + ] + class GeneralL2Lasso(Optimizer): r"""The Minimizer class solves the following equation, @@ -144,31 +277,12 @@ class GeneralL2Lasso(Optimizer): alpha: float = 1e-3 lambda1: float = 1e-6 hyperparameters: dict = None - inverse_dimension: List[ - Union[cp.Dimension, cp.LinearDimension, cp.MonotonicDimension] - ] = [] - max_iterations: int = 10000 - tolerance: float = 1e-5 - positive: bool = True - regularizer: Literal["smooth lasso", "sparse ridge fusion"] = "smooth lasso" - method: Literal["multi-task", "gradient_decent", "lars"] = "gradient_decent" - f: Union[cp.CSDM, np.ndarray] = None - n_iter: int = None - _scale: float = PrivateAttr(1.0) - _estimator: Any = PrivateAttr() def __init__(self, **kwargs): super().__init__(**kwargs) if self.hyperparameters is None: self.hyperparameters = {"lambda": self.lambda1, "alpha": self.alpha} - class Config: - arbitrary_types_allowed = True - - @property - def f_shape(self): - return tuple([item.count for item in self.inverse_dimension])[::-1] - def fit(self, K, s): r"""Fit the model using the coordinate descent method from scikit-learn. @@ -183,52 +297,16 @@ def fit(self, K, s): :math:`{\bf s}`, as a :math:`m \times m_\text{count}` matrix. """ - s_ = _get_proper_data(s) - self._scale = s_.real.max() - s_ = s_ / self._scale + Ks, ss, (half, s_) = self._fit_setup(K, s, self.hyperparameters["alpha"], False) - # prod = np.asarray(self.f_shape).prod() - # if K.shape[1] != prod: - # raise ValueError( - # "The product of the shape, `f_shape`, must be equal to the length of " - # f"the axis 1 of kernel, K, {K.shape[1]} != {prod}." - # ) - - half = False - if "half" in self.inverse_dimension[0].application: - half = self.inverse_dimension[0].application["half"] - - Ks, ss = _get_augmented_data( - K=K, - s=s_, - alpha=s_.size * self.hyperparameters["alpha"], - regularizer=self.regularizer, - f_shape=self.f_shape, - half=half, - ) - - _estimator = self._get_minimizer(self.hyperparameters["lambda"]) + _estimator = self.get_optimizer(self.hyperparameters["lambda"]) _estimator.fit(Ks, ss) f = _estimator.coef_.copy() - if half: - index = self.inverse_dimension[0].application["index"] - f_new = np.zeros((s_.shape[1], np.prod(self.f_shape)), dtype=float) - f_new[:, index] = f - f = f_new - if s_.shape[1] > 1: - f.shape = (s_.shape[1],) + self.f_shape - f[:, :, 0] /= 2.0 - f[:, 0, :] /= 2.0 - else: - f.shape = self.f_shape - f[:, 0] /= 2.0 - f[0, :] /= 2.0 - - f *= self._scale + f = self.shape_solution(f, s_, half) if isinstance(s, cp.CSDM): - f = self.pack_as_csdm(f, s) + f = self.pack_solution_as_csdm(f, s) # f = cp.as_csdm(f) # if len(s.x) > 1: @@ -240,86 +318,13 @@ def fit(self, K, s): self.f = f self.n_iter = _estimator.n_iter_ - def pack_as_csdm(self, f, s): - f = cp.as_csdm(f) - - if len(s.x) > 1: - f.x[2] = s.x[1] - f.x[1] = self.inverse_dimension[1] - f.x[0] = self.inverse_dimension[0] - return f - - def predict(self, K): - r"""Predict the signal using the linear model. - - Args - ---- - - K: ndarray - A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of shape - (m, n). - - Return - ------ - ndarray - A numpy array of shape (m, m_count) with the predicted values - """ - return self._estimator.predict(K) * self._scale - - def residuals(self, K, s): - r"""Return the residual as the difference the data and the prediced data(fit), - following - - .. math:: - \text{residuals} = {\bf s - Kf^*} - - where :math:`{\bf f^*}` is the optimum solution. - - Args - ---- - K: ndarray. - A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of shape - (m, n). - s: ndarray ot CSDM object. - A csdm object or a :math:`m \times m_\text{count}` signal matrix, - :math:`{\bf s}`. - - Return - ------ - ndarray or CSDM object. - If `s` is a csdm object, returns a csdm object with the residuals. If `s` - is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. - csdm object - """ - s_ = s.y[0].components[0].T if isinstance(s, cp.CSDM) else s - predict = np.squeeze(self._estimator.predict(K)) * self._scale - residue = s_ - predict - - if not isinstance(s, cp.CSDM): - return residue - - residue = cp.as_csdm(residue.T) - residue._dimensions = s._dimensions - return residue - - def score(self, K, s, sample_weights=None): - """ - The coefficient of determination, :math:`R^2`, of the prediction. - For more information, read scikit-learn documentation. - """ - return self._estimator.score(K, s / self._scale, sample_weights) - class GeneralL2LassoLS(GeneralL2Lasso): - folds: int = 10 - sigma: float = 0.0 - randomize: bool = False - times: int = 2 verbose: bool = False n_jobs: int = -1 _path: Any = PrivateAttr() - def fit(self, K, s, **kwargs): + def fit(self, K, s, n_repeats=1, **kwargs): 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. @@ -329,75 +334,74 @@ def fit(self, K, s, **kwargs): shape (m, n). s: A :math:`m \times m_\text{count}` signal matrix, :math:`{\bf s}` as a csdm object or a numpy array or shape (m, m_count). + n_repeats: int. Repeat the cross-validation by adding noise to the input + data and evaluate the mean cross-validation error. """ - Ks, ss, cv_indexes = self._pre_fit_cv(K, s, self.hyperparameters["alpha"]) + Ks, ss, cv_indexes = self._fit_setup( + K, s, self.hyperparameters["alpha"], get_cv_indexes=True + ) + l1 = self.get_optimizer(self.hyperparameters["lambda"]) + + ks0 = K.shape[0] + new_indexes = self._get_bootstrap_cv_indexes(cv_indexes, n_repeats, ks0) self._path = [] sigma_sq = self.sigma ** 2 - ks0 = K.shape[0] - l1 = self._get_minimizer(self.hyperparameters["lambda"]) - - def fnc(x0): - self._path.append(x0) - Ks[ks0:] *= np.sqrt(x0[0] / self.hyperparameters["alpha"]) - alpha = self.hyperparameters["alpha"] = x0[0] - l1.alpha = x0[1] / 2.0 - mse = -cv(l1, Ks, ss, cv_indexes, alpha=alpha, n_jobs=-1) - mse -= sigma_sq - # mse *= 1000 - return np.sqrt(2 * mse) - # + np.mean( - # (l1.coef_[1:, :] - l1.coef_[:-1, :]) ** 2 - # ) - # for i in range(self.folds): - # # print(i) - # fit_ = l1.fit(Ks_[cv_indexes[i][0]], ss[cv_indexes[i][0]]) - # y_predict = fit_.predict(Ks_[cv_indexes[i][1]]) - # avg += np.sum((y_predict - ss[cv_indexes[i][1]]) ** 2) - # return avg + def fnc(x_): + x = [10 ** x_[0], 10 ** x_[1]] + self._path.append(x) + Ks[ks0:] *= np.sqrt(x[0] / self.hyperparameters["alpha"]) + l1.alpha = x[1] / 2.0 - x0 = [self.hyperparameters["alpha"], self.hyperparameters["lambda"]] + mse = np.sum( + [cv(l1, Ks, ss, indexes, n_jobs=self.n_jobs) for indexes in new_indexes] + ) + mse /= -n_repeats + mse -= sigma_sq + return np.sqrt(2 * mse) - bounds = [(1e-15, 1e-15), (1e-4, 1e-4)] - # method="Powell" - res = least_squares(fnc, x0, bounds=bounds, verbose=2, **kwargs) + x0 = np.log10([self.hyperparameters["alpha"], self.hyperparameters["lambda"]]) + + kwargs["xtol"] = 1e-4 if "xtol" not in kwargs else kwargs["xtol"] + res = least_squares( + fnc, + x0, + # jac="3-point", + # loss="cauchy", + # method="dogbox", + bounds=[(-15, -15), (-4, -4)], + # diff_step=[5e-5, 5e-5], + verbose=2, + **kwargs, + ) - self.hyperparameters["alpha"] = res.x[0] - self.hyperparameters["lambda"] = res.x[1] + self.hyperparameters["alpha"] = 10 ** res.x[0] + self.hyperparameters["lambda"] = 10 ** res.x[1] - print(res.message, res.status) + print(res.message) + print(res.status) # Calculate the solution using the complete data at the optimized lambda and # alpha values super().fit(K, s) + @property + def path(self): + alpha_ = [item[0] for item in self._path] + lambda_ = [item[1] for item in self._path] + return {"alpha": alpha_, "lambda": lambda_} -class GeneralL2LassoCV(Optimizer): +class GeneralL2LassoCV(Optimizer): alphas: Union[List[float], np.ndarray] = None lambdas: Union[List[float], np.ndarray] = None hyperparameters: dict = None - inverse_dimension: List[ - Union[cp.Dimension, cp.LinearDimension, cp.MonotonicDimension] - ] = [] - max_iterations: int = 10000 - tolerance: float = 1e-5 - positive: bool = True - regularizer: Literal["smooth lasso", "sparse ridge fusion"] = "smooth lasso" - method: Literal["multi-task", "gradient_decent", "lars"] = "gradient_decent" - folds: int = 10 - sigma: float = 0.0 - randomize: bool = False - times: int = 2 verbose: bool = False n_jobs: int = -1 + f_std: Union[cp.CSDM, np.ndarray] = None - f: Union[cp.CSDM, np.ndarray] = None - n_iter: int = None - _scale: float = PrivateAttr(1.0) - _opt: Any = PrivateAttr() _path: Any = PrivateAttr() _cv_alphas: Any = PrivateAttr() _cv_lambdas: Any = PrivateAttr() @@ -406,10 +410,6 @@ class GeneralL2LassoCV(Optimizer): class Config: arbitrary_types_allowed = True - @property - def f_shape(self): - return tuple([item.count for item in self.inverse_dimension])[::-1] - def __init__(self, **kwargs): super().__init__(**kwargs) self._cv_alphas = ( @@ -426,78 +426,93 @@ def __init__(self, **kwargs): if self.hyperparameters is None: self.hyperparameters = {"lambda": None, "alpha": None} - def fit(self, K, s, scoring="neg_mean_squared_error"): - 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. + def get_optimizer_jobs(self, Ks, ss, cv_indexes, start_index, n_repeats=1): + """Generate a grid of n_alpha x n_lambda optimizers (jobs) over the given + grid of hyperparameters. Args: - K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). - s: A :math:`m \times m_\text{count}` signal matrix, :math:`{\bf s}` as a - csdm object or a numpy array or shape (m, m_count). + np.ndarray Ks: The augmented kernel + cp.CSDM ss: The augmented dataset + (list of tuple) cv_indexes: (test, train) indexes for cross-validation + int start_index: The index in the augmented kernel (Ks), where the + smoothening matrix starts. + int n_repeats: Repeat the cross-validation by adding noise to the input + data and evaluate the mean cross-validation error. """ + # create an l1 over all lambda values + l1 = self.get_optimizer(self._cv_lambdas[0]) - Ks, ss, cv_indexes = self._pre_fit_cv(K, s, self._cv_alphas[0]) - - self._cv_map = np.zeros((self._cv_alphas.size, self._cv_lambdas.size)) + l1_array = [] + for lambda_ in self._cv_lambdas: + l1_array += [deepcopy(l1)] + l1_array[-1].alpha = lambda_ / 2.0 + # l2 optimizer per l1 over alpha values alpha_ratio = np.ones(self._cv_alphas.size) if self._cv_alphas.size != 1 and self._cv_alphas[0] != 0: alpha_ratio[1:] = np.sqrt(self._cv_alphas[1:] / self._cv_alphas[:-1]) - start_index = K.shape[0] - - l1 = self._get_minimizer(self._cv_lambdas[0]) - l1_array = [] - - for lambda_ in self._cv_lambdas: - l1_array.append(deepcopy(l1)) - l1_array[-1].alpha = lambda_ / 2.0 + def alpha_lambda_grid_jobs(new_indexes): + """Create optimization jobs over alpha-lambda grid + Args: + ss: Augmented dataset. + """ + jobs_ = [] + for alpha_ratio_ in alpha_ratio: + if alpha_ratio_ != 0: + Ks[start_index:] *= alpha_ratio_ + Ks_copy = deepcopy(Ks) + jobs_ += [ + delayed(cv)(l1_, Ks_copy, ss, new_indexes) for l1_ in l1_array + ] + return jobs_ + + new_indexes = self._get_bootstrap_cv_indexes(cv_indexes, n_repeats, start_index) + return [ + item for indexes in new_indexes for item in alpha_lambda_grid_jobs(indexes) + ] + + def fit(self, K, s, scoring="neg_mean_squared_error", n_repeats=1): + 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. - j = 0 - for alpha_ratio_ in alpha_ratio: - if alpha_ratio_ != 0: - Ks[start_index:] *= alpha_ratio_ - jobs = ( - delayed(cv)(l1_, Ks, ss, cv_indexes, alpha=alpha_) - for l1_, alpha_ in zip(l1_array, self._cv_alphas) - ) - self._cv_map[j] = Parallel( - n_jobs=self.n_jobs, - verbose=self.verbose, - **{ - "backend": { - "threads": "threading", - "processes": "multiprocessing", - None: None, - }["threads"] - }, - )(jobs) - j += 1 - - # _cv_map contains negated mean square errors, therefore multiply by -1. - self._cv_map *= -1 - # subtract the variance. - if scoring == "neg_mean_squared_error": - self._cv_map -= self.sigma ** 2 - index = np.where(self._cv_map < 0) - self._cv_map[index] = np.nan + Args: + K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of + shape (m, n). + s: A :math:`m \times m_\text{count}` signal matrix, :math:`{\bf s}` as a + csdm object or a numpy array or shape (m, m_count). + n_repeats: int. Repeat the cross-validation by adding noise to the input + data and evaluate the mean cross-validation error. + """ + # get augmented kernel, dataset, and cv-indexes. + Ks, ss, cv_indexes = self._fit_setup(K, s, self._cv_alphas[0], True) + + # optimizers as jobs over the hyperparameter grid. + jobs = self.get_optimizer_jobs(Ks, ss, cv_indexes, K.shape[0], n_repeats) + + # start parallel processes for running jobs + processes = Parallel( + n_jobs=self.n_jobs, + verbose=self.verbose, + **{ + "backend": { + "threads": "threading", + "processes": "multiprocessing", + None: None, + }["threads"] + }, + ) + result = processes(jobs) # run jobs - # After subtracting the variance, any negative values in the cv grid is a - # result of fitting noise. Take the absolute value of cv to avoid such - # models. - # self._cv_map = np.abs(self._cv_map) + cv_map = self.get_cv_map(result, n_repeats) # get cross-validation error metric - # The argmin of the minimum value is the selected model as it has the least - # prediction error. - index = np.unravel_index(self._cv_map.argmin(), self._cv_map.shape) - self.hyperparameters["alpha"] = self._cv_alphas[index[0]] - self.hyperparameters["lambda"] = self._cv_lambdas[index[1]] + # find the minima of the curve to get the optimum hyperparameters. + self.find_optimum_hyperparameters(cv_map, scoring) # Calculate the solution using the complete data at the optimized lambda and - # alpha values - self._opt = GeneralL2Lasso( + # alpha values. + opt = GeneralL2Lasso( alpha=self.hyperparameters["alpha"], lambda1=self.hyperparameters["lambda"], max_iterations=self.max_iterations, @@ -507,62 +522,65 @@ def fit(self, K, s, scoring="neg_mean_squared_error"): inverse_dimension=self.inverse_dimension, method=self.method, ) - self._opt.fit(K, s) - self.f = self._opt.f + opt.fit(K, s) + self._estimator = opt._estimator + self.f = opt.f + self.n_iter = opt.n_iter - # convert _cv_map to csdm - self._cv_map = cp.as_csdm(np.squeeze(self._cv_map.T)) - if len(self._cv_alphas) != 1: - d0 = cp.as_dimension(-np.log10(self._cv_alphas), label="-log(α)") - self._cv_map.x[0] = d0 + def get_cv_map(self, cv_map, n_repeats=1, scoring="neg_mean_squared_error"): + # cv_map = np.asarray([item["test_score"].mean() for item in result]) + cv_map = np.asarray(cv_map) + cv_map.shape = (n_repeats, self._cv_alphas.size, self._cv_lambdas.size) - if len(self._cv_lambdas) == 1: - return + # cv_map contains negated mean square errors, therefore multiply by -1. + cv_map *= -1 - d1 = cp.as_dimension(-np.log10(self._cv_lambdas), label="-log(λ)") - if len(self._cv_alphas) != 1: - self._cv_map.x[1] = d1 - else: - self._cv_map.x[0] = d1 + # subtract the variance. + if scoring == "neg_mean_squared_error": + cv_map -= self.sigma ** 2 - def predict(self, K): - r"""Predict the signal using the linear model. + # After subtracting the variance, any negative values in the cv grid is a + # result of noise. Substitute negative values for np.nan to avoid such + # models. + index = np.where(cv_map < 0) + cv_map[index] = np.nan - Args: - K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). + return cv_map.mean(axis=0) - Return: - A numpy array of shape (m, m_count) with the predicted values. - """ - return self._opt.predict(K) + def find_optimum_hyperparameters(self, cv_map, scoring="neg_mean_squared_error"): + """Find the global minimal of the cross-validation error metric and update + the hyperparameter values. - def residuals(self, K, s): - r"""Return the residual as the difference the data and the prediced data(fit), - following + Args: + np.nd.array cv_map: cross-calidation error metric array. + str scoring: Scoring criteria for error metric. + """ - .. math:: - \text{residuals} = {\bf s - Kf^*} + # The argmin of the minimum value is the selected model as it has the least + # prediction error. + cv_shape = (self._cv_alphas.size, self._cv_lambdas.size) + index = np.unravel_index(cv_map.argmin(), cv_shape) + self.hyperparameters["alpha"] = self._cv_alphas[index[0]] + self.hyperparameters["lambda"] = self._cv_lambdas[index[1]] + self.pack_cv_map_as_csdm(cv_map) - where :math:`{\bf f^*}` is the optimum solution. + def pack_cv_map_as_csdm(self, cv_map): + """Convert the cross-validation error metric as a CSDM object. Args: - K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). - s: A csdm object or a :math:`m \times m_\text{count}` signal matrix, - :math:`{\bf s}`. - Return: - If `s` is a csdm object, returns a csdm object with the residuals. If `s` - is a numpy array, return a :math:`m \times m_\text{count}` residue matrix. + np.nd.array cv_map: cross-calidation error metric array. """ - return self._opt.residuals(K, s) + cv_map = cp.as_csdm(cv_map.T) + cv_map.x[0] = cp.as_dimension(-np.log10(self._cv_alphas), label="-log(α)") + cv_map.x[1] = cp.as_dimension(-np.log10(self._cv_lambdas), label="-log(λ)") - def score(self, K, s, sample_weights=None): - """ - Return the coefficient of determination, :math:`R^2`, of the prediction. - For more information, read scikit-learn documentation. - """ - return self._opt.score(K, s, sample_weights) + if self._cv_alphas.size == 1: + cv_map = cv_map[0] + + if self._cv_lambdas.size == 1: + cv_map = cv_map[:, 0] + + self._cv_map = cv_map @property def cross_validation_curve(self): @@ -573,7 +591,7 @@ def cross_validation_curve(self): return self._cv_map -def cv(l1, X, y, cv, alpha=0, n_jobs=1): +def cv(l1, X, y, cv_indexes, n_jobs=1): """Return the cross-validation score as negative of mean square error.""" if isinstance(l1, (Lasso, MultiTaskLasso)): fit_params = {"check_input": False} @@ -585,11 +603,11 @@ def cv(l1, X, y, cv, alpha=0, n_jobs=1): X=X, y=y, scoring="neg_mean_absolute_error", - cv=cv, + cv=cv_indexes, fit_params=fit_params, n_jobs=n_jobs, verbose=0, - return_estimator=True, + # return_estimator=True, ) # estimators = cv_score["estimator"] # print(estimators) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..d26b730 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,6 @@ +[project] +requires-python = ">=3.6" + +[tool.pylint.messages_control] +extension-pkg-whitelist = ["pydantic", "mrsimulator"] +disable = "R0201" From 830b3d9b42139679d144e010c5b6246ed8a16734 Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Mon, 26 Apr 2021 12:50:43 -0400 Subject: [PATCH 10/11] fixed examples. --- .../auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb | 2 +- .../synthetic/plot_1D_2_sideband_bimodal.ipynb | 2 +- .../auto_examples/synthetic/plot_1D_3_MAF.ipynb | 10 +++++----- examples/sideband/plot_2D_KMg0.5O4SiO2.py | 4 ++-- examples/synthetic/plot_1D_3_MAF.py | 5 ++--- 5 files changed, 11 insertions(+), 12 deletions(-) diff --git a/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb b/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb index f73fc0c..c73819e 100644 --- a/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb +++ b/docs/notebooks/auto_examples/sideband/plot_2D_KMg0.5O4SiO2.ipynb @@ -220,7 +220,7 @@ }, "outputs": [], "source": [ - "# setup the pre-defined range of alpha and lambda values\nlambdas = 10 ** (-4.6 - 1 * (np.arange(5) / 4))\nalphas = 10 ** (-4.8 - 1.5 * (np.arange(5) / 4))\n\n# setup the smooth lasso cross-validation class\ns_lasso = SmoothLassoCV(\n alphas=alphas, # A numpy array of alpha values.\n lambdas=lambdas, # A numpy array of lambda values.\n sigma=0.00070, # The standard deviation of noise from the MAT dataset.\n folds=10, # The number of folds in n-folds cross-validation.\n inverse_dimension=inverse_dimensions, # previously defined inverse dimensions.\n verbose=1, # If non-zero, prints the progress as the computation proceeds.\n max_iterations=20000, # The maximum number of allowed interations.\n)\n\n# run fit using the compressed kernel and compressed data.\ns_lasso.fit(compressed_K, compressed_s)" + "# setup the pre-defined range of alpha and lambda values\nlambdas = 10 ** (-5.4 - 1 * (np.arange(5) / 4))\nalphas = 10 ** (-4.5 - 1.5 * (np.arange(5) / 4))\n\n# setup the smooth lasso cross-validation class\ns_lasso = SmoothLassoCV(\n alphas=alphas, # A numpy array of alpha values.\n lambdas=lambdas, # A numpy array of lambda values.\n sigma=0.00070, # The standard deviation of noise from the MAT dataset.\n folds=10, # The number of folds in n-folds cross-validation.\n inverse_dimension=inverse_dimensions, # previously defined inverse dimensions.\n verbose=1, # If non-zero, prints the progress as the computation proceeds.\n max_iterations=20000, # The maximum number of allowed interations.\n)\n\n# run fit using the compressed kernel and compressed data.\ns_lasso.fit(compressed_K, compressed_s)" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb index a9a8bd4..98e480e 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_2_sideband_bimodal.ipynb @@ -231,7 +231,7 @@ }, "outputs": [], "source": [ - "lambdas = 10 ** (-5 - 1 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n # max_iterations=50000,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s, n_repeats=2)" + "lambdas = 10 ** (-5 - 1 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" ] }, { diff --git a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb index 565b438..4045687 100644 --- a/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb +++ b/docs/notebooks/auto_examples/synthetic/plot_1D_3_MAF.ipynb @@ -33,7 +33,7 @@ }, "outputs": [], "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import TSVDCompression, SmoothLasso, SmoothLassoCV\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.set_xlim(0, 100)\n ax.set_ylim(0, 100)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom pylab import rcParams\n\nfrom mrinversion.kernel.nmr import ShieldingPALineshape\nfrom mrinversion.linear_model import SmoothLasso, SmoothLassoCV, TSVDCompression\nfrom mrinversion.utils import get_polar_grids\n\n# Setup for the matplotlib figures\nrcParams[\"figure.figsize\"] = 4.5, 3.5\nrcParams[\"font.size\"] = 9\n\n\n# function for 2D x-y plot.\ndef plot2D(ax, csdm_object, title=\"\"):\n # convert the dimension coordinates of the csdm_object from Hz to pmm.\n csdm_object.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n csdm_object.dimensions[1].to(\"ppm\", \"nmr_frequency_ratio\")\n\n levels = (np.arange(9) + 1) / 10\n ax.contourf(csdm_object, cmap=\"gist_ncar\", levels=levels)\n ax.set_xlim(0, 100)\n ax.set_ylim(0, 100)\n ax.set_title(title)\n get_polar_grids(ax)\n ax.set_aspect(\"equal\")" ] }, { @@ -249,7 +249,7 @@ }, "outputs": [], "source": [ - "lambdas = 10 ** (-5.2 - 1.5 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 1.5 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n max_iterations=50000,\n)\n\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" + "lambdas = 10 ** (-4.8 - 1.5 * (np.arange(6) / 5))\nalphas = 10 ** (-4 - 2 * (np.arange(6) / 5))\n\n# set up cross validation smooth lasso method\ns_lasso_cv = SmoothLassoCV(\n alphas=alphas,\n lambdas=lambdas,\n inverse_dimension=inverse_dimension,\n sigma=0.005,\n folds=10,\n)\n# run the fit using the compressed kernel and compressed signal.\ns_lasso_cv.fit(compressed_K, compressed_s)" ] }, { @@ -285,14 +285,14 @@ }, "outputs": [], "source": [ - "cv_curve = s_lasso_cv.cross_validation_curve\n\n# plot of the cross-validation curve\nplt.figure(figsize=(5, 3.5))\nax = plt.subplot(projection=\"csdm\")\nax.contour(np.log10(s_lasso_cv.cross_validation_curve), levels=25)\nax.scatter(\n -np.log10(s_lasso_cv.hyperparameters[\"alpha\"]),\n -np.log10(s_lasso_cv.hyperparameters[\"lambda\"]),\n marker=\"x\",\n color=\"k\",\n label=\"Grid evaluation\",\n)\nplt.legend()\nplt.tight_layout(pad=0.5)\nplt.show()" + "cv_curve = s_lasso_cv.cross_validation_curve\n\n# plot of the cross-validation curve\nplt.figure(figsize=(5, 3.5))\nax = plt.subplot(projection=\"csdm\")\nax.contour(np.log10(s_lasso_cv.cross_validation_curve), levels=25)\nax.scatter(\n -np.log10(s_lasso_cv.hyperparameters[\"alpha\"]),\n -np.log10(s_lasso_cv.hyperparameters[\"lambda\"]),\n marker=\"x\",\n color=\"k\",\n)\nplt.tight_layout(pad=0.5)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### The optimum solution\n\nThe :attr:`~mrinversion.linear_model.SmoothLassoLS.f` attribute of the instance holds\nthe solution.\n\n" + "### The optimum solution\n\nThe :attr:`~mrinversion.linear_model.SmoothLassoCV.f` attribute of the instance holds\nthe solution.\n\n" ] }, { @@ -303,7 +303,7 @@ }, "outputs": [], "source": [ - "f_sol = s_lasso.f" + "f_sol = s_lasso_cv.f" ] }, { diff --git a/examples/sideband/plot_2D_KMg0.5O4SiO2.py b/examples/sideband/plot_2D_KMg0.5O4SiO2.py index d9e9f4c..8ef494c 100644 --- a/examples/sideband/plot_2D_KMg0.5O4SiO2.py +++ b/examples/sideband/plot_2D_KMg0.5O4SiO2.py @@ -195,8 +195,8 @@ def plot2D(csdm_object, **kwargs): # increase the grid resolution for your problem if desired. # setup the pre-defined range of alpha and lambda values -lambdas = 10 ** (-4.6 - 1 * (np.arange(5) / 4)) -alphas = 10 ** (-4.8 - 1.5 * (np.arange(5) / 4)) +lambdas = 10 ** (-5.4 - 1 * (np.arange(5) / 4)) +alphas = 10 ** (-4.5 - 1.5 * (np.arange(5) / 4)) # setup the smooth lasso cross-validation class s_lasso = SmoothLassoCV( diff --git a/examples/synthetic/plot_1D_3_MAF.py b/examples/synthetic/plot_1D_3_MAF.py index 46f9671..a5b852d 100644 --- a/examples/synthetic/plot_1D_3_MAF.py +++ b/examples/synthetic/plot_1D_3_MAF.py @@ -216,8 +216,8 @@ def plot2D(ax, csdm_object, title=""): # hyperparameters. # The following code generates a range of :math:`\lambda` and :math:`\alpha` values # that are uniformly sampled on the log scale. -lambdas = 10 ** (-5.2 - 1.5 * (np.arange(6) / 5)) -alphas = 10 ** (-4 - 1.5 * (np.arange(6) / 5)) +lambdas = 10 ** (-4.8 - 1.5 * (np.arange(6) / 5)) +alphas = 10 ** (-4 - 2 * (np.arange(6) / 5)) # set up cross validation smooth lasso method s_lasso_cv = SmoothLassoCV( @@ -259,7 +259,6 @@ def plot2D(ax, csdm_object, title=""): marker="x", color="k", ) -plt.legend() plt.tight_layout(pad=0.5) plt.show() From 3ad7be63850ed7d5b5ef401422dd243bc8e1f56b Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava Date: Fri, 29 Oct 2021 18:02:02 -0400 Subject: [PATCH 11/11] fixed T2 kernel generation. --- mrinversion/kernel/base.py | 17 ++++++++++++++++- mrinversion/kernel/relaxation.py | 21 +++++++++------------ mrinversion/linear_model/_base_l1l2.py | 11 ++++++++--- 3 files changed, 33 insertions(+), 16 deletions(-) diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 541d121..6eb2804 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -15,12 +15,27 @@ class BaseModel: """Base kernel class.""" - def __init__(self, kernel_dimension, inverse_kernel_dimension, n_dir, n_inv): + def __init__( + self, + kernel_dimension: list, + inverse_kernel_dimension: list, + n_dir: int, + n_inv: int, + ): + """Initialize the class + + Args: + kernel_dimension: A list of kernel CSDM dimensions. + inverse_kernel_dimension: A list of inverse kernel CSDM dimensions. + n_dir: Number of kernel dimensions. + n_inv: Number of inverse dimensions. + """ kernel = self.__class__.__name__ message = ( f"Exactly {n_inv} inverse dimension(s) is/are required for the " f"{kernel} kernel." ) + if isinstance(inverse_kernel_dimension, list): if len(inverse_kernel_dimension) != n_inv: raise ValueError(message) diff --git a/mrinversion/kernel/relaxation.py b/mrinversion/kernel/relaxation.py index 1604d86..1a28894 100644 --- a/mrinversion/kernel/relaxation.py +++ b/mrinversion/kernel/relaxation.py @@ -6,8 +6,7 @@ class T2(BaseModel): - r""" - A class for simulating the kernel of T2 decaying functions, + r"""A class for simulating the kernel of T2 decaying functions, .. math:: y = \exp(-x/x_\text{inv}). @@ -24,29 +23,27 @@ def __init__(self, kernel_dimension, inverse_dimension): # inverse_dimension = inverse_dimension * cp.ScalarQuantity("s") super().__init__(kernel_dimension, inverse_dimension, 1, 1) - def kernel(self): - """ - Return the kernel of T2 decaying functions. + def kernel(self, supersampling=1): + """Return the kernel of T2 decaying functions. Returns: A numpy array. """ x = self.kernel_dimension.coordinates - x_inverse = self.inverse_kernel_dimension.coordinates + # x_inverse = self.inverse_kernel_dimension.coordinates # x = self.kernel_dimension.coordinates.to("s").value # self.inverse_kernel_dimension = ( # self.inverse_kernel_dimension / cp.ScalarQuantity("s") # ) - # x_inverse = _supersampled_coordinates( - # self.inverse_kernel_dimension, supersampling=supersampling - # ) + x_inverse = _supersampled_coordinates( + self.inverse_kernel_dimension, supersampling=supersampling + ) amp = np.exp(np.tensordot(-(1 / x_inverse), x, 0)) - return self._averaged_kernel(amp, 1) + return self._averaged_kernel(amp, supersampling) class T1(BaseModel): - r""" - A class for simulating the kernel of T1 recovery functions, + r"""A class for simulating the kernel of T1 recovery functions, .. math:: y = 1 - \exp(-x/x_\text{inv}). diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 0898021..6fcf3e6 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -170,13 +170,18 @@ def pack_solution_as_csdm(self, f, s): """ f = cp.as_csdm(f) + for i, item in enumerate(self.inverse_dimension): + f.x[i] = item + if len(s.x) > 1: - f.x[2] = s.x[1] - f.x[1] = self.inverse_dimension[1] - f.x[0] = self.inverse_dimension[0] + f.x[i + 1] = s.x[1] return f def shape_solution(self, f, s_, half=False): + if len(self.f_shape) == 1: + f *= self._scale + return f + if half: index = self.inverse_dimension[0].application["index"] f_new = np.zeros((s_.shape[1], np.prod(self.f_shape)), dtype=float)