diff --git a/.flake8 b/.flake8 index 0909542..1b6de64 100644 --- a/.flake8 +++ b/.flake8 @@ -1,9 +1,9 @@ -[flake8] -ignore = E203 -exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py -filename = *.pyx, *.px*, *py -max-line-length = 88 -max-complexity = 10 -select = C,E,F,W,N8 -count = True -statistics = True +[flake8] +ignore = E203 +exclude = .eggs, *.egg,build, mrinversion/kernel/__init__.py, mrinversion/utils.py, examples/* +filename = *.pyx, *.px*, *py +max-line-length = 88 +max-complexity = 10 +select = C,E,F,W,N8 +count = True +statistics = True diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index 51ffa22..56f24e4 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -1,58 +1,58 @@ -# For most projects, this workflow file will not need changing; you simply need -# to commit it to your repository. -# -# You may wish to alter this file to override the set of languages analyzed, -# or to provide custom queries or build logic. -name: "CodeQL" - -on: - # At the end of every day - schedule: - - cron: '0 0 * * *' - -jobs: - analyze: - name: Analyze - runs-on: ubuntu-latest - - strategy: - fail-fast: false - matrix: - # Override automatic language detection by changing the below list - # Supported options are ['csharp', 'cpp', 'go', 'java', 'javascript', 'python'] - language: ['python'] - # Learn more... - # https://docs.github.com/en/github/finding-security-vulnerabilities-and-errors-in-your-code/configuring-code-scanning#overriding-automatic-language-detection - - steps: - - name: Checkout repository - uses: actions/checkout@v6 - - # Initializes the CodeQL tools for scanning. - - name: Initialize CodeQL - uses: github/codeql-action/init@v4 - with: - languages: ${{ matrix.language }} - # If you wish to specify custom queries, you can do so here or in a config file. - # By default, queries listed here will override any specified in a config file. - # Prefix the list here with "+" to use these queries and those in the config file. - # queries: ./path/to/local/query, your-org/your-repo/queries@main - - # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). - # If this step fails, then you should remove it and run the build manually (see below) - - name: Autobuild - uses: github/codeql-action/autobuild@v4 - - # ℹī¸ Command-line programs to run using the OS shell. - # 📚 https://git.io/JvXDl - - # ✏ī¸ If the Autobuild fails above, remove it and uncomment the following three lines - # and modify them (or add more) to build your code if your project - # uses a compiled language - - #- run: | - # make bootstrap - # make release - - - name: Perform CodeQL Analysis - uses: github/codeql-action/analyze@v4 +# For most projects, this workflow file will not need changing; you simply need +# to commit it to your repository. +# +# You may wish to alter this file to override the set of languages analyzed, +# or to provide custom queries or build logic. +name: "CodeQL" + +on: + # At the end of every day + schedule: + - cron: '0 0 * * *' + +jobs: + analyze: + name: Analyze + runs-on: ubuntu-latest + + strategy: + fail-fast: false + matrix: + # Override automatic language detection by changing the below list + # Supported options are ['csharp', 'cpp', 'go', 'java', 'javascript', 'python'] + language: ['python'] + # Learn more... + # https://docs.github.com/en/github/finding-security-vulnerabilities-and-errors-in-your-code/configuring-code-scanning#overriding-automatic-language-detection + + steps: + - name: Checkout repository + uses: actions/checkout@v6 + + # Initializes the CodeQL tools for scanning. + - name: Initialize CodeQL + uses: github/codeql-action/init@v4 + with: + languages: ${{ matrix.language }} + # If you wish to specify custom queries, you can do so here or in a config file. + # By default, queries listed here will override any specified in a config file. + # Prefix the list here with "+" to use these queries and those in the config file. + # queries: ./path/to/local/query, your-org/your-repo/queries@main + + # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). + # If this step fails, then you should remove it and run the build manually (see below) + - name: Autobuild + uses: github/codeql-action/autobuild@v4 + + # ℹī¸ Command-line programs to run using the OS shell. + # 📚 https://git.io/JvXDl + + # ✏ī¸ If the Autobuild fails above, remove it and uncomment the following three lines + # and modify them (or add more) to build your code if your project + # uses a compiled language + + #- run: | + # make bootstrap + # make release + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v4 diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 277fb2d..0df86f3 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -1,66 +1,66 @@ -# This workflow will install Python dependencies, run tests and lint with a variety of Python versions -# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions - -name: CI - -on: - push: - branches: [master] - pull_request: - branches: [master] - -jobs: - code_lint: - runs-on: ubuntu-latest - strategy: - matrix: - python-version: ["3.11"] - - steps: - - name: Checkout - uses: actions/checkout@v6 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v6 - with: - python-version: ${{ matrix.python-version }} - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install flake8 - - name: Lint with flake8 - run: | - # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=B,C,E,F,W,T,N8 --ignore=E203 --max-line-length=88 --show-source --statistics --exclude="examples/* docs/* build/* utils.py" - # 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/* docs/* build/* utils.py" - - test_os: - needs: [code_lint] - runs-on: "ubuntu-latest" - strategy: - matrix: - python-version: ["3.10", "3.11", "3.12", "3.13"] - - steps: - - name: Checkout - uses: actions/checkout@v6 - - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v6 - with: - python-version: ${{ matrix.python-version }} - - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install setuptools - pip install -r requirements-dev.txt - - - name: Build and install package from source - run: python setup.py develop - - - name: Test with pytest - run: pytest --cov=./ --cov-report=xml - - - name: Upload coverage - uses: codecov/codecov-action@v5.5.1 +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: CI + +on: + push: + branches: [master] + pull_request: + branches: [master] + +jobs: + code_lint: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.11"] + + steps: + - name: Checkout + uses: actions/checkout@v6 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install flake8 + - name: Lint with flake8 + run: | + # stop the build if there are Python syntax errors or undefined names + flake8 . --count --select=B,C,E,F,W,T,N8 --ignore=E203 --max-line-length=88 --show-source --statistics --exclude="examples/* docs/* build/* utils.py" + # 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/* docs/* build/* utils.py" + + test_os: + needs: [code_lint] + runs-on: "ubuntu-latest" + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13"] + + steps: + - name: Checkout + uses: actions/checkout@v6 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install setuptools + pip install -r requirements-dev.txt + + - name: Build and install package from source + run: python setup.py develop + + - name: Test with pytest + run: pytest --cov=./ --cov-report=xml + + - name: Upload coverage + uses: codecov/codecov-action@v5.5.1 diff --git a/CHANGELOG b/CHANGELOG index ea9f43a..d353611 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,30 +1,30 @@ -v0.3.1 ------- - -What's new! -''''''''''' - -- Simplified plot_3d function -- Replace fortran code with numba jit python version. - -v0.3.0 ------- - -What's new! -''''''''''' - -- Added T1, T2 relaxation inversion kernel. - -v0.2.0 ------- - -What's new! -''''''''''' - -- Added :func:`~mrinversion.utils.to_Haeberlen_grid` function to convert the 3D :math:`\rho(\delta_\text{iso}, x, y)` - distribution to :math:`\rho(\delta_\text{iso}, \zeta_\sigma, \eta_\sigma)` distribution. - -Changes -''''''' - -- Update code to comply with updated csdmpy library. +v0.3.1 +------ + +What's new! +''''''''''' + +- Simplified plot_3d function +- Replace fortran code with numba jit python version. + +v0.3.0 +------ + +What's new! +''''''''''' + +- Added T1, T2 relaxation inversion kernel. + +v0.2.0 +------ + +What's new! +''''''''''' + +- Added :func:`~mrinversion.utils.to_Haeberlen_grid` function to convert the 3D :math:`\rho(\delta_\text{iso}, x, y)` + distribution to :math:`\rho(\delta_\text{iso}, \zeta_\sigma, \eta_\sigma)` distribution. + +Changes +''''''' + +- Update code to comply with updated csdmpy library. diff --git a/docs/conf.py b/docs/conf.py index be2e5f0..f52a298 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,362 +1,362 @@ -# Configuration file for the Sphinx documentation builder. -# -# This file does only contain a selection of the most common options. For a -# full list see the documentation: -# http://www.sphinx-doc.org/en/master/config -# -- Path setup -------------------------------------------------------------- -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# -import os -import sys -import warnings - -from sphinx_gallery.sorting import ExplicitOrder -from sphinx_gallery.sorting import FileNameSortKey - -sys.path.insert(0, os.path.abspath("../..")) - -# -- Project information ----------------------------------------------------- -project = "mrinversion" -copyright = "2020, Deepansh J. Srivastava" -author = "Deepansh J. Srivastava" - -path = os.path.split(__file__)[0] -# get version number from the file -with open(os.path.join(path, "../mrinversion/__init__.py")) as f: - for line in f.readlines(): - if "__version__" in line: - before_keyword, keyword, after_keyword = line.partition("=") - __version__ = after_keyword.strip()[1:-1] - -# The short X.Y version -version = __version__ -# The full version, including alpha/beta/rc tags -release = __version__ - - -# -- General configuration --------------------------------------------------- - -# If your documentation needs a minimal Sphinx version, state it here. -needs_sphinx = "4.0" - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. -extensions = [ - "sphinx.ext.autodoc", - "matplotlib.sphinxext.plot_directive", - "sphinx.ext.doctest", - "sphinx.ext.mathjax", - "sphinx.ext.autosummary", - "sphinx.ext.viewcode", - "sphinx.ext.napoleon", - "sphinx_copybutton", - "sphinxjp.themes.basicstrap", - "sphinx_gallery.gen_gallery", - "sphinx.ext.intersphinx", - "sphinx_tabs.tabs", -] - -autosummary_generate = True - -# ---------------------------------------------------------------------------- # -# Plot directive config # -# ---------------------------------------------------------------------------- # -plot_html_show_source_link = False -plot_rcparams = { - "font.size": 10, - "font.weight": "light", - "font.family": "sans-serif", - "font.sans-serif": "Helvetica", -} - -# ---------------------------------------------------------------------------- # -# Sphinx Gallery config # -# ---------------------------------------------------------------------------- # - -# filter sphinx matplotlib warning -warnings.filterwarnings( - "ignore", - category=UserWarning, - message="Matplotlib is currently using agg, which is a" - " non-GUI backend, so cannot show the figure.", -) - -# numfig config -numfig = True -numfig_secnum_depth = 1 -numfig_format = {"figure": "Figure %s", "table": "Table %s", "code-block": "Listing %s"} - -# math -math_number_all = True - -# sphinx gallery config -sphinx_gallery_conf = { - "examples_dirs": ["../examples"], # path to example scripts - "remove_config_comments": True, - "gallery_dirs": ["galley_examples"], # path to gallery generated output - "within_subsection_order": FileNameSortKey, - "subsection_order": ExplicitOrder( - [ - "../examples/synthetic", - "../examples/sideband", - "../examples/MAF", - "../examples/relaxation", - ] - ), - "reference_url": { - # The module you locally document uses None - "mrinversion": None, - }, - "first_notebook_cell": ( - "# This cell is added by sphinx-gallery\n\n" - "%matplotlib inline\n\n" - "import mrinversion\n" - "print(f'You are using mrinversion v{mrinversion.__version__}')" - ), - # "binder": { - # # Required keys - # "org": "DeepanshS", - # "repo": "mrinversion", - # "branch": "master", - # "binderhub_url": "https://mybinder.org", - # "dependencies": "../requirements.txt", - # # Optional keys - # "filepath_prefix": "docs/_build/html", - # "notebooks_dir": "../../notebooks", - # "use_jupyter_lab": True, - # }, -} - -intersphinx_mapping = { - "matplotlib": ("https://matplotlib.org", None), - "numpy": ("https://numpy.org/doc/stable/", None), - "csdmpy": ("https://csdmpy.readthedocs.io/en/latest/", None), - "astropy": ("https://docs.astropy.org/en/stable/", None), -} - -copybutton_prompt_text = ">>> |\\$ |\\[\\d*\\]: |\\.\\.\\.: |[.][.][.] " -copybutton_prompt_is_regexp = True - -# ---------------------------------------------------------------------------- # - -# Add any paths that contain templates here, relative to this directory. -templates_path = ["_templates"] - -# The suffix(es) of source filenames. -# You can specify multiple suffix as a list of string: -source_suffix = ".rst" - -# The master toc-tree document. -master_doc = "index" - -# autodoc mock modules -autodoc_mock_imports = [] - -# The language for content autogenerated by Sphinx. Refer to documentation -# for a list of supported languages. -# -# This is also used if you do content translation via gettext catalogs. -# Usually you set "language" from the command line for these cases. -language = "en" - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -# This pattern also affects html_static_path and html_extra_path. -exclude_patterns = ["_build", "**.ipynb_checkpoints", "Thumbs.db", ".DS_Store"] - -# The name of the Pygments (syntax highlighting) style to use. -pygments_style = None - - -# -- Options for HTML output ------------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# - -# Some html_theme options are 'alabaster', 'bootstrap', 'sphinx_rtd_theme', -# 'classic', 'basicstrap' -html_theme = "basicstrap" - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. -# -html_theme_options = { - # Set the lang attribute of the html tag. Defaults to 'en' - "lang": "en", - # Disable showing the sidebar. Defaults to 'false' - "nosidebar": False, - # Show header searchbox. Defaults to false. works only "nosidebar=True", - "header_searchbox": False, - # Put the sidebar on the right side. Defaults to false. - "rightsidebar": False, - # Set the width of the sidebar. Defaults to 3 - "sidebar_span": 3, - # Fix navbar to top of screen. Defaults to true - "nav_fixed_top": True, - # Fix the width of the sidebar. Defaults to false - "nav_fixed": True, - # Set the width of the sidebar. Defaults to '900px' - "nav_width": "300px", - # Fix the width of the content area. Defaults to false - "content_fixed": False, - # Set the width of the content area. Defaults to '900px' - "content_width": "900px", - # Fix the width of the row. Defaults to false - "row_fixed": False, - # Disable the responsive design. Defaults to false - "noresponsive": False, - # Disable the responsive footer relbar. Defaults to false - "noresponsiverelbar": False, - # Disable flat design. Defaults to false. - # Works only "bootstrap_version = 3" - "noflatdesign": False, - # Enable Google Web Font. Defaults to false - "googlewebfont": False, - # Set the URL of Google Web Font's CSS. - # Defaults to 'http://fonts.googleapis.com/css?family=Text+Me+One' - # "googlewebfont_url": "http://fonts.googleapis.com/css?family=Roboto", # NOQA - # Set the Style of Google Web Font's CSS. - # Defaults to "font-family: 'Text Me One', sans-serif;" - # "googlewebfont_style": u"font-family: 'Roboto' Regular;", # font-size: 1.5em", - # Set 'navbar-inverse' attribute to header navbar. Defaults to false. - "header_inverse": True, - # Set 'navbar-inverse' attribute to relbar navbar. Defaults to false. - "relbar_inverse": True, - # Enable inner theme by Bootswatch. Defaults to false - "inner_theme": False, - # Set the name of inner theme. Defaults to 'bootswatch-simplex' - # "inner_theme_name": "bootswatch-simplex", - # Select Twitter bootstrap version 2 or 3. Defaults to '3' - "bootstrap_version": "3", - # Show "theme preview" button in header navbar. Defaults to false. - "theme_preview": False, - # Set the Size of Heading text. Defaults to None - # "h1_size": "3.0em", - # "h2_size": "2.6em", - # "h3_size": "2.2em", - # "h4_size": "1.8em", - # "h5_size": "1.9em", - # "h6_size": "1.1em", -} - - -# Theme options -html_logo = "_static/mrinversion.png" -html_style = "style.css" -html_title = f"mrinversion:doc v{__version__}" -html_last_updated_fmt = "" -# html_logo = "mrinversion" -html_sidebars = { - "**": ["searchbox.html", "globaltoc.html"], - "using/windows": ["searchbox.html", "windowssidebar.html"], -} - - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ["_static"] - - -# -- Options for HTMLHelp output --------------------------------------------- - -# Output file base name for HTML help builder. -htmlhelp_basename = "Mrinversion doc" - -# -- Options for LaTeX output ------------------------------------------------ -latex_engine = "xelatex" -# latex_logo = "_static/csdmpy.png" -latex_show_pagerefs = True - -latex_elements = { - # The paper size ('letterpaper' or 'a4paper'). - # - "papersize": "letterpaper", - # The font size ('10pt', '11pt' or '12pt'). - # - "pointsize": "9pt", - "fontenc": r"\usepackage[utf8]{inputenc}", - "geometry": r"\usepackage[vmargin=2.5cm, hmargin=2cm]{geometry}", - # "fncychap": r"\usepackage[Rejne]{fncychap}", - # Additional stuff for the LaTeX preamble. - "preamble": r""" - \usepackage[T1]{fontenc} - \usepackage{amsfonts, amsmath, amssymb, mathbbol} - \usepackage{graphicx} - \usepackage{setspace} - \singlespacing - - \usepackage{fancyhdr} - \pagestyle{fancy} - \fancyhf{} - \fancyhead[L]{ - \ifthenelse{\isodd{\value{page}}}{ \small \nouppercase{\leftmark} }{} - } - \fancyhead[R]{ - \ifthenelse{\isodd{\value{page}}}{}{ \small \nouppercase{\rightmark} } - } - \fancyfoot[CO, CE]{\thepage} - """, - # Latex figure (float) alignment - # - "figure_align": "htbp", -} - -# Grouping the document tree into LaTeX files. List of tuples -# (source start file, target name, title, -# author, documentclass [howto, manual, or own class]). -latex_documents = [ - (master_doc, "mrinversion.tex", "mrinversion Documentation", author, "manual") -] - - -# -- Options for manual page output ------------------------------------------ - -# One entry per manual page. List of tuples -# (source start file, name, description, authors, manual section). -man_pages = [(master_doc, "mrinversion", "mrinversion Documentation", [author], 1)] - - -# -- Options for Texinfo output ---------------------------------------------- - -# Grouping the document tree into Texinfo files. List of tuples -# (source start file, target name, title, author, -# dir menu entry, description, category) -texinfo_documents = [ - ( - master_doc, - "Mrinversion", - "Mrinversion Documentation", - author, - "Mrinversion", - "Statistical learning of tensor distribution from NMR anisotropic spectra", - "Miscellaneous", - ) -] - - -# -- Options for Epub output ------------------------------------------------- - -# Bibliographic Dublin Core info. -epub_title = project - -# The unique identifier of the text. This can be a ISBN number -# or the project homepage. -# -# epub_identifier = '' - -# A unique identification for the text. -# -# epub_uid = '' - -# A list of files that should not be packed into the epub file. -epub_exclude_files = ["search.html", "_static/style.css"] - - -def setup(app): - app.add_css_file("style.css") +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config +# -- Path setup -------------------------------------------------------------- +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +import warnings + +from sphinx_gallery.sorting import ExplicitOrder +from sphinx_gallery.sorting import FileNameSortKey + +sys.path.insert(0, os.path.abspath("../..")) + +# -- Project information ----------------------------------------------------- +project = "mrinversion" +copyright = "2020, Deepansh J. Srivastava" +author = "Deepansh J. Srivastava" + +path = os.path.split(__file__)[0] +# get version number from the file +with open(os.path.join(path, "../mrinversion/__init__.py")) as f: + for line in f.readlines(): + if "__version__" in line: + before_keyword, keyword, after_keyword = line.partition("=") + __version__ = after_keyword.strip()[1:-1] + +# The short X.Y version +version = __version__ +# The full version, including alpha/beta/rc tags +release = __version__ + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +needs_sphinx = "4.0" + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + "sphinx.ext.autodoc", + "matplotlib.sphinxext.plot_directive", + "sphinx.ext.doctest", + "sphinx.ext.mathjax", + "sphinx.ext.autosummary", + "sphinx.ext.viewcode", + "sphinx.ext.napoleon", + "sphinx_copybutton", + "sphinxjp.themes.basicstrap", + "sphinx_gallery.gen_gallery", + "sphinx.ext.intersphinx", + "sphinx_tabs.tabs", +] + +autosummary_generate = True + +# ---------------------------------------------------------------------------- # +# Plot directive config # +# ---------------------------------------------------------------------------- # +plot_html_show_source_link = False +plot_rcparams = { + "font.size": 10, + "font.weight": "light", + "font.family": "sans-serif", + "font.sans-serif": "Helvetica", +} + +# ---------------------------------------------------------------------------- # +# Sphinx Gallery config # +# ---------------------------------------------------------------------------- # + +# filter sphinx matplotlib warning +warnings.filterwarnings( + "ignore", + category=UserWarning, + message="Matplotlib is currently using agg, which is a" + " non-GUI backend, so cannot show the figure.", +) + +# numfig config +numfig = True +numfig_secnum_depth = 1 +numfig_format = {"figure": "Figure %s", "table": "Table %s", "code-block": "Listing %s"} + +# math +math_number_all = True + +# sphinx gallery config +sphinx_gallery_conf = { + "examples_dirs": ["../examples"], # path to example scripts + "remove_config_comments": True, + "gallery_dirs": ["galley_examples"], # path to gallery generated output + "within_subsection_order": FileNameSortKey, + "subsection_order": ExplicitOrder( + [ + "../examples/synthetic", + "../examples/sideband", + "../examples/MAF", + "../examples/relaxation", + ] + ), + "reference_url": { + # The module you locally document uses None + "mrinversion": None, + }, + "first_notebook_cell": ( + "# This cell is added by sphinx-gallery\n\n" + "%matplotlib inline\n\n" + "import mrinversion\n" + "print(f'You are using mrinversion v{mrinversion.__version__}')" + ), + # "binder": { + # # Required keys + # "org": "DeepanshS", + # "repo": "mrinversion", + # "branch": "master", + # "binderhub_url": "https://mybinder.org", + # "dependencies": "../requirements.txt", + # # Optional keys + # "filepath_prefix": "docs/_build/html", + # "notebooks_dir": "../../notebooks", + # "use_jupyter_lab": True, + # }, +} + +intersphinx_mapping = { + "matplotlib": ("https://matplotlib.org", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "csdmpy": ("https://csdmpy.readthedocs.io/en/latest/", None), + "astropy": ("https://docs.astropy.org/en/stable/", None), +} + +copybutton_prompt_text = ">>> |\\$ |\\[\\d*\\]: |\\.\\.\\.: |[.][.][.] " +copybutton_prompt_is_regexp = True + +# ---------------------------------------------------------------------------- # + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +source_suffix = ".rst" + +# The master toc-tree document. +master_doc = "index" + +# autodoc mock modules +autodoc_mock_imports = [] + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = "en" + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ["_build", "**.ipynb_checkpoints", "Thumbs.db", ".DS_Store"] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = None + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# + +# Some html_theme options are 'alabaster', 'bootstrap', 'sphinx_rtd_theme', +# 'classic', 'basicstrap' +html_theme = "basicstrap" + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +html_theme_options = { + # Set the lang attribute of the html tag. Defaults to 'en' + "lang": "en", + # Disable showing the sidebar. Defaults to 'false' + "nosidebar": False, + # Show header searchbox. Defaults to false. works only "nosidebar=True", + "header_searchbox": False, + # Put the sidebar on the right side. Defaults to false. + "rightsidebar": False, + # Set the width of the sidebar. Defaults to 3 + "sidebar_span": 3, + # Fix navbar to top of screen. Defaults to true + "nav_fixed_top": True, + # Fix the width of the sidebar. Defaults to false + "nav_fixed": True, + # Set the width of the sidebar. Defaults to '900px' + "nav_width": "300px", + # Fix the width of the content area. Defaults to false + "content_fixed": False, + # Set the width of the content area. Defaults to '900px' + "content_width": "900px", + # Fix the width of the row. Defaults to false + "row_fixed": False, + # Disable the responsive design. Defaults to false + "noresponsive": False, + # Disable the responsive footer relbar. Defaults to false + "noresponsiverelbar": False, + # Disable flat design. Defaults to false. + # Works only "bootstrap_version = 3" + "noflatdesign": False, + # Enable Google Web Font. Defaults to false + "googlewebfont": False, + # Set the URL of Google Web Font's CSS. + # Defaults to 'http://fonts.googleapis.com/css?family=Text+Me+One' + # "googlewebfont_url": "http://fonts.googleapis.com/css?family=Roboto", # NOQA + # Set the Style of Google Web Font's CSS. + # Defaults to "font-family: 'Text Me One', sans-serif;" + # "googlewebfont_style": u"font-family: 'Roboto' Regular;", # font-size: 1.5em", + # Set 'navbar-inverse' attribute to header navbar. Defaults to false. + "header_inverse": True, + # Set 'navbar-inverse' attribute to relbar navbar. Defaults to false. + "relbar_inverse": True, + # Enable inner theme by Bootswatch. Defaults to false + "inner_theme": False, + # Set the name of inner theme. Defaults to 'bootswatch-simplex' + # "inner_theme_name": "bootswatch-simplex", + # Select Twitter bootstrap version 2 or 3. Defaults to '3' + "bootstrap_version": "3", + # Show "theme preview" button in header navbar. Defaults to false. + "theme_preview": False, + # Set the Size of Heading text. Defaults to None + # "h1_size": "3.0em", + # "h2_size": "2.6em", + # "h3_size": "2.2em", + # "h4_size": "1.8em", + # "h5_size": "1.9em", + # "h6_size": "1.1em", +} + + +# Theme options +html_logo = "_static/mrinversion.png" +html_style = "style.css" +html_title = f"mrinversion: doc v{__version__}" +html_last_updated_fmt = "" +# html_logo = "mrinversion" +html_sidebars = { + "**": ["searchbox.html", "globaltoc.html"], + "using/windows": ["searchbox.html", "windowssidebar.html"], +} + + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = "Mrinversion doc" + +# -- Options for LaTeX output ------------------------------------------------ +latex_engine = "xelatex" +# latex_logo = "_static/csdmpy.png" +latex_show_pagerefs = True + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + "papersize": "letterpaper", + # The font size ('10pt', '11pt' or '12pt'). + # + "pointsize": "9pt", + "fontenc": r"\usepackage[utf8]{inputenc}", + "geometry": r"\usepackage[vmargin=2.5cm, hmargin=2cm]{geometry}", + # "fncychap": r"\usepackage[Rejne]{fncychap}", + # Additional stuff for the LaTeX preamble. + "preamble": r""" + \usepackage[T1]{fontenc} + \usepackage{amsfonts, amsmath, amssymb, mathbbol} + \usepackage{graphicx} + \usepackage{setspace} + \singlespacing + + \usepackage{fancyhdr} + \pagestyle{fancy} + \fancyhf{} + \fancyhead[L]{ + \ifthenelse{\isodd{\value{page}}}{ \small \nouppercase{\leftmark} }{} + } + \fancyhead[R]{ + \ifthenelse{\isodd{\value{page}}}{}{ \small \nouppercase{\rightmark} } + } + \fancyfoot[CO, CE]{\thepage} + """, + # Latex figure (float) alignment + # + "figure_align": "htbp", +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, "mrinversion.tex", "mrinversion Documentation", author, "manual") +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [(master_doc, "mrinversion", "mrinversion Documentation", [author], 1)] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + ( + master_doc, + "Mrinversion", + "Mrinversion Documentation", + author, + "Mrinversion", + "Statistical learning of tensor distribution from NMR anisotropic spectra", + "Miscellaneous", + ) +] + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ["search.html", "_static/style.css"] + + +def setup(app): + app.add_css_file("style.css") diff --git a/mrinversion/__init__.py b/mrinversion/__init__.py index 260c070..6285450 100644 --- a/mrinversion/__init__.py +++ b/mrinversion/__init__.py @@ -1 +1 @@ -__version__ = "0.3.1" +__version__ = "0.3.1" diff --git a/mrinversion/kernel/__init__.py b/mrinversion/kernel/__init__.py index d6c5c68..0a6d164 100644 --- a/mrinversion/kernel/__init__.py +++ b/mrinversion/kernel/__init__.py @@ -1,3 +1,3 @@ -"""Initialization of the kernel module.""" -from mrinversion.kernel.relaxation import T1 # NOQA -from mrinversion.kernel.relaxation import T2 # NOQA +"""Initialization of the kernel module.""" +from mrinversion.kernel.relaxation import T1 # NOQA +from mrinversion.kernel.relaxation import T2 # NOQA diff --git a/mrinversion/kernel/base.py b/mrinversion/kernel/base.py index 681d95c..ad16f9c 100644 --- a/mrinversion/kernel/base.py +++ b/mrinversion/kernel/base.py @@ -1,205 +1,219 @@ -"""Base classes for kernel generation.""" -from copy import deepcopy - -import csdmpy as cp -import numpy as np -from mrsimulator.method.lib import BlochDecaySpectrum - -from .utils import _x_y_to_zeta_eta_distribution - -__dimension_list__ = (cp.Dimension, cp.LinearDimension, cp.MonotonicDimension) -__dimension_name__ = ("Dimension", "LinearDimension", "MonotonicDimension") - - -class BaseModel: - """Base kernel class.""" - - def __init__(self, kernel_dimension, inverse_kernel_dimension, n_dir, n_inv): - 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) - if n_inv == 1: - inverse_kernel_dimension = inverse_kernel_dimension[0] - else: - if n_inv != 1: - raise ValueError(message) - - inverse_kernel_dimension = _check_csdm_dimension( - inverse_kernel_dimension, "inverse_dimension" - ) - - message = ( - f"Exactly {n_dir} direct dimension(s) is/are required for the " - f"{kernel} kernel." - ) - if isinstance(kernel_dimension, list): - if len(kernel_dimension) != n_dir: - raise ValueError(message) - if n_dir == 1: - kernel_dimension = kernel_dimension[0] - else: - if n_dir != 1: - raise ValueError(message) - - kernel_dimension = _check_csdm_dimension(kernel_dimension, "kernel_dimension") - - self.kernel_dimension = kernel_dimension - self.inverse_kernel_dimension = inverse_kernel_dimension - self.inverse_dimension = self.inverse_kernel_dimension - - def _averaged_kernel(self, amp, supersampling, xy_grid=True): - """Return the kernel by averaging over the supersampled grid cells.""" - shape = () - inverse_kernel_dimension = self.inverse_kernel_dimension - if not isinstance(self.inverse_kernel_dimension, list): - inverse_kernel_dimension = [self.inverse_kernel_dimension] - - for item in inverse_kernel_dimension[::-1]: - shape += (item.count, supersampling) - shape += (self.kernel_dimension.count,) - - K = amp.reshape(shape) - - inv_len = len(inverse_kernel_dimension) - axes = tuple(2 * i + 1 for i in range(inv_len)) - K = K.mean(axis=axes) - - if xy_grid: - section = [*[0 for _ in range(inv_len)], slice(None, None, None)] - K /= K[tuple(section)].sum() - - section = [slice(None, None, None) for _ in range(inv_len + 1)] - for i, item in enumerate(inverse_kernel_dimension): - if item.coordinates[0].value == 0: - section_ = deepcopy(section) - section_[i] = 0 - K[tuple(section_)] /= 2.0 - - inv_size = np.asarray([item.count for item in inverse_kernel_dimension]).prod() - K = K.reshape(inv_size, self.kernel_dimension.count).T - return K - - -class LineShape(BaseModel): - """Base line-shape kernel generation class.""" - - def __init__( - self, - kernel_dimension, - inverse_kernel_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency=None, - number_of_sidebands=None, - ): - super().__init__(kernel_dimension, inverse_kernel_dimension, 1, 2) - - kernel = self.__class__.__name__ - dim_types = ["frequency", "dimensionless"] - _check_dimension_type(self.kernel_dimension, "anisotropic", dim_types, kernel) - _check_dimension_type( - self.inverse_kernel_dimension, "inverse", dim_types, kernel - ) - - dim = self.kernel_dimension - - temp_method = BlochDecaySpectrum.parse_dict_with_units( - {"channels": [channel], "magnetic_flux_density": magnetic_flux_density} - ) - # larmor frequency from method. - B0 = temp_method.spectral_dimensions[0].events[0].magnetic_flux_density # in T - gamma = temp_method.channels[0].gyromagnetic_ratio # in MHz/T - self.larmor_frequency = -gamma * B0 # in MHz - - spectral_width = dim.increment * dim.count - reference_offset = dim.coordinates_offset - if dim.complex_fft is False: - reference_offset = dim.coordinates_offset + spectral_width / 2.0 - - if dim.increment.unit.physical_type == "dimensionless": - lf = abs(self.larmor_frequency) - val = spectral_width.to("ppm") - spectral_width = f"{np.round(val.value * lf, decimals=8)} Hz" - val = reference_offset.to("ppm") - reference_offset = f"{np.round(val.value * lf, decimals=8)} Hz" - - spectral_dimensions = [ - dict( - count=dim.count, - reference_offset=str(reference_offset), - spectral_width=str(spectral_width), - ) - ] - - if rotor_frequency is None: - if dim.increment.unit.physical_type == "dimensionless": - rotor_frequency = f"{dim.increment.to('ppm').value * lf} Hz" - else: - rotor_frequency = str(dim.increment) - - self.method_args = { - "channels": [channel], - "magnetic_flux_density": magnetic_flux_density, - "rotor_angle": rotor_angle, - "rotor_frequency": rotor_frequency, - "spectral_dimensions": spectral_dimensions, - } - - self.number_of_sidebands = number_of_sidebands - if number_of_sidebands is None: - self.number_of_sidebands = dim.count - - def _get_zeta_eta(self, supersampling): - """Return zeta and eta coordinates over x-y grid""" - - zeta, eta = _x_y_to_zeta_eta_distribution( - self.inverse_kernel_dimension, supersampling - ) - return zeta, eta - - -def _check_csdm_dimension(dimensions, dimension_id): - if not isinstance(dimensions, (list, *__dimension_list__)): - raise ValueError( - f"The value of the `{dimension_id}` attribute must be one of " - f"`{__dimension_name__}` objects." - ) - - # copy the list - # dimensions = deepcopy(dimensions) - if isinstance(dimensions, __dimension_list__): - return dimensions - - for i, item in enumerate(dimensions): - if not isinstance(item, __dimension_list__): - raise ValueError( - f"The element at index {i} of the `{dimension_id}` list must be an " - f"instance of one of the {__dimension_name__} classes." - ) - # dimensions[i] = item - return dimensions - - -def _check_dimension_type(dimensions, direction, dimension_quantity, kernel_type): - if not isinstance(dimensions, list): - if dimensions.quantity_name not in dimension_quantity: - raise ValueError( - f"A {direction} dimension with quantity name `{dimension_quantity}` " - f"is required for the `{kernel_type}` kernel, instead got " - f"`{dimensions.quantity_name}` as the quantity name for the dimension." - ) - return - for i, item in enumerate(dimensions): - if item.quantity_name not in dimension_quantity: - raise ValueError( - f"A {direction} dimension with quantity name `{dimension_quantity}` " - f"is required for the `{kernel_type}` kernel, instead got " - f"`{item.quantity_name}` as the quantity name for the dimension at " - f"index {i}." - ) +"""Base classes for kernel generation.""" +from copy import deepcopy + +import csdmpy as cp +import numpy as np +from mrsimulator.method.lib import BlochDecaySpectrum + +from .utils import _x_y_to_zeta_eta_distribution + +__dimension_list__ = (cp.Dimension, cp.LinearDimension, cp.MonotonicDimension) +__dimension_name__ = ("Dimension", "LinearDimension", "MonotonicDimension") + + +class BaseModel: + """Base kernel class.""" + + def __init__(self, kernel_dimension, inverse_kernel_dimension, n_dir, n_inv): + 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) + if n_inv == 1: + inverse_kernel_dimension = inverse_kernel_dimension[0] + else: + if n_inv != 1: + raise ValueError(message) + + inverse_kernel_dimension = _check_csdm_dimension( + inverse_kernel_dimension, "inverse_dimension" + ) + + message = ( + f"Exactly {n_dir} direct dimension(s) is/are required for the " + f"{kernel} kernel." + ) + if isinstance(kernel_dimension, list): + if len(kernel_dimension) != n_dir: + raise ValueError(message) + if n_dir == 1: + kernel_dimension = kernel_dimension[0] + else: + if n_dir != 1: + raise ValueError(message) + + kernel_dimension = _check_csdm_dimension(kernel_dimension, "kernel_dimension") + + self.kernel_dimension = kernel_dimension + self.inverse_kernel_dimension = inverse_kernel_dimension + self.inverse_dimension = self.inverse_kernel_dimension + + def _averaged_kernel(self, amp, supersampling, xy_grid=True, mask_kernel=False): + """Return the kernel by averaging over the supersampled grid cells.""" + shape = () + inverse_kernel_dimension = self.inverse_kernel_dimension + if not isinstance(self.inverse_kernel_dimension, list): + inverse_kernel_dimension = [self.inverse_kernel_dimension] + + for item in inverse_kernel_dimension[::-1]: + shape += (item.count, supersampling) + shape += (self.kernel_dimension.count,) + # print(f'shape: {shape}') + K = amp.reshape(shape) + + inv_len = len(inverse_kernel_dimension) + axes = tuple(2 * i + 1 for i in range(inv_len)) + K = K.mean(axis=axes) + # print(f'K shape: {K.shape}') + + if xy_grid: + section = [*[0 for _ in range(inv_len)], slice(None, None, None)] + # print(f'section: {section}') + # print(f'K[tuple(section)]: {K[tuple(section)]}') + K /= K[tuple(section)].sum() + + section = [slice(None, None, None) for _ in range(inv_len + 1)] + for i, item in enumerate(inverse_kernel_dimension): + if item.coordinates[0].value == 0: + section_ = deepcopy(section) + section_[i] = 0 + K[tuple(section_)] /= 2.0 + # print(f'K-shape: {K.shape}') + # if mask_kernel: + # mask = np.zeros(K.shape) + # for i in range(len(K)): + # mask[i][len(K[0])-i:] += 1 + # K = ma.masked_array(K, mask=mask) + + inv_size = np.asarray([item.count for item in inverse_kernel_dimension]).prod() + K = K.reshape(inv_size, self.kernel_dimension.count).T + return K + + +class LineShape(BaseModel): + """Base line-shape kernel generation class.""" + + def __init__( + self, + kernel_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency=None, + number_of_sidebands=None, + ): + super().__init__(kernel_dimension, inverse_kernel_dimension, 1, 2) + + kernel = self.__class__.__name__ + dim_types = ["frequency", "dimensionless"] + _check_dimension_type(self.kernel_dimension, "anisotropic", dim_types, kernel) + _check_dimension_type( + self.inverse_kernel_dimension, "inverse", dim_types, kernel + ) + + dim = self.kernel_dimension + + temp_method = BlochDecaySpectrum.parse_dict_with_units( + {"channels": [channel], "magnetic_flux_density": magnetic_flux_density} + ) + # larmor frequency from method. + B0 = temp_method.spectral_dimensions[0].events[0].magnetic_flux_density # in T + gamma = temp_method.channels[0].gyromagnetic_ratio # in MHz/T + self.larmor_frequency = -gamma * B0 # in MHz + + spectral_width = dim.increment * dim.count + reference_offset = dim.coordinates_offset + if dim.complex_fft is False: + reference_offset = dim.coordinates_offset + spectral_width / 2.0 + + if dim.increment.unit.physical_type == "dimensionless": + lf = abs(self.larmor_frequency) + val = spectral_width.to("ppm") + spectral_width = f"{np.round(val.value * lf, decimals=8)} Hz" + val = reference_offset.to("ppm") + reference_offset = f"{np.round(val.value * lf, decimals=8)} Hz" + + spectral_dimensions = [ + dict( + count=dim.count, + reference_offset=str(reference_offset), + spectral_width=str(spectral_width), + ) + ] + + if rotor_frequency is None: + if dim.increment.unit.physical_type == "dimensionless": + rotor_frequency = f"{dim.increment.to('ppm').value * lf} Hz" + else: + rotor_frequency = str(dim.increment) + + self.method_args = { + "channels": [channel], + "magnetic_flux_density": magnetic_flux_density, + "rotor_angle": rotor_angle, + "rotor_frequency": rotor_frequency, + "spectral_dimensions": spectral_dimensions, + } + + self.number_of_sidebands = number_of_sidebands + if number_of_sidebands is None: + self.number_of_sidebands = dim.count + + def _get_zeta_eta(self, supersampling, eta_bound=1, calc_pos=False): + """Return zeta and eta coordinates over x-y grid""" + if eta_bound == 1 and not calc_pos: + zeta, eta = _x_y_to_zeta_eta_distribution( + self.inverse_kernel_dimension, supersampling, eta_bound + ) + return zeta, eta + else: + zeta, eta, abundances = _x_y_to_zeta_eta_distribution( + self.inverse_kernel_dimension, supersampling, eta_bound, calc_pos + ) + return zeta, eta, abundances + + +def _check_csdm_dimension(dimensions, dimension_id): + if not isinstance(dimensions, (list, *__dimension_list__)): + raise ValueError( + f"The value of the `{dimension_id}` attribute must be one of " + f"`{__dimension_name__}` objects." + ) + + # copy the list + # dimensions = deepcopy(dimensions) + if isinstance(dimensions, __dimension_list__): + return dimensions + + for i, item in enumerate(dimensions): + if not isinstance(item, __dimension_list__): + raise ValueError( + f"The element at index {i} of the `{dimension_id}` list must be an " + f"instance of one of the {__dimension_name__} classes." + ) + # dimensions[i] = item + return dimensions + + +def _check_dimension_type(dimensions, direction, dimension_quantity, kernel_type): + if not isinstance(dimensions, list): + if dimensions.quantity_name not in dimension_quantity: + raise ValueError( + f"A {direction} dimension with quantity name `{dimension_quantity}` " + f"is required for the `{kernel_type}` kernel, instead got " + f"`{dimensions.quantity_name}` as the quantity name for the dimension." + ) + return + for i, item in enumerate(dimensions): + if item.quantity_name not in dimension_quantity: + raise ValueError( + f"A {direction} dimension with quantity name `{dimension_quantity}` " + f"is required for the `{kernel_type}` kernel, instead got " + f"`{item.quantity_name}` as the quantity name for the dimension at " + f"index {i}." + ) diff --git a/mrinversion/kernel/csa_aniso.py b/mrinversion/kernel/csa_aniso.py index a27dcf8..4d47ebd 100644 --- a/mrinversion/kernel/csa_aniso.py +++ b/mrinversion/kernel/csa_aniso.py @@ -1,223 +1,181 @@ -from copy import deepcopy - -from mrsimulator import Simulator -from mrsimulator import SpinSystem -from mrsimulator.method.lib import BlochDecaySpectrum - -from mrinversion.kernel.base import LineShape - - -class ShieldingPALineshape(LineShape): - """ - A generalized class for simulating the pure anisotropic NMR nuclear shielding - line-shape kernel. - - Args: - anisotropic_dimension: A Dimension object, or an equivalent dictionary - object. This dimension must represent the pure anisotropic - dimension. - inverse_dimension: A list of two Dimension objects, or equivalent - dictionary objects representing the `x`-`y` coordinate grid. - channel: The channel is an isotope symbol of the nuclei given as the atomic - number followed by the atomic symbol, for example, `1H`, `13C`, and - `29Si`. This nucleus must correspond to the recorded frequency - resonances. - magnetic_flux_density: The magnetic flux density of the external static - magnetic field. The default value is 9.4 T. - rotor_angle: The angle of the sample holder (rotor) relative to the - direction of the external magnetic field. The default value is - 54.735 deg (magic angle). - rotor_frequency: The effective sample spin rate. Depending on the NMR - sequence, this value may be less than the physical sample rotation - frequency. The default is 14 kHz. - number_of_sidebands: The number of sidebands to simulate along the - anisotropic dimension. The default value is 1. - """ - - def __init__( - self, - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density="9.4 T", - rotor_angle="54.735 deg", - rotor_frequency="14 kHz", - number_of_sidebands=1, - ): - super().__init__( - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density, - rotor_angle, - rotor_frequency, - number_of_sidebands, - ) - - def kernel(self, supersampling=1): - """ - Return the NMR nuclear shielding anisotropic line-shape kernel. - - Args: - supersampling: An integer. Each cell is supersampled by the factor - `supersampling` along every dimension. - Returns: - A numpy array containing the line-shape kernel. - """ - args_ = deepcopy(self.method_args) - method = BlochDecaySpectrum.parse_dict_with_units(args_) - isotope = args_["channels"][0] - zeta, eta = self._get_zeta_eta(supersampling) - - x_csdm = self.inverse_kernel_dimension[0] - if x_csdm.coordinates.unit.physical_type == "frequency": - # convert zeta to ppm if given in frequency units. - zeta /= self.larmor_frequency # zeta in ppm - - for dim_i in self.inverse_kernel_dimension: - if dim_i.origin_offset.value == 0: - dim_i.origin_offset = f"{abs(self.larmor_frequency)} MHz" - - spin_systems = [ - SpinSystem( - sites=[dict(isotope=isotope, shielding_symmetric=dict(zeta=z, eta=e))] - ) - for z, e in zip(zeta, eta) - ] - - sim = Simulator() - sim.config.number_of_sidebands = self.number_of_sidebands - sim.config.decompose_spectrum = "spin_system" - - sim.spin_systems = spin_systems - sim.methods = [method] - sim.run(pack_as_csdm=False) - - amp = sim.methods[0].simulation.real - return self._averaged_kernel(amp, supersampling) - - -class MAF(ShieldingPALineshape): - r"""A specialized class for simulating the pure anisotropic NMR nuclear shielding - line-shape kernel resulting from the 2D MAF spectra. - - Args: - anisotropic_dimension: A Dimension object, or an equivalent dictionary - object. This dimension must represent the pure anisotropic - dimension. - inverse_dimension: A list of two Dimension objects, or equivalent - dictionary objects representing the `x`-`y` coordinate grid. - channel: The isotope symbol of the nuclei given as the atomic number - followed by the atomic symbol, for example, `1H`, `13C`, and - `29Si`. This nucleus must correspond to the recorded frequency - resonances. - magnetic_flux_density: The magnetic flux density of the external static - magnetic field. The default value is 9.4 T. - - **Assumptions:** - The simulated line-shapes correspond to an infinite speed spectrum spinning at - :math:`90^\circ`. - """ - - def __init__( - self, - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density="9.4 T", - ): - - super().__init__( - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density, - "90 deg", - "200 GHz", - 1, - ) - - -class SpinningSidebands(ShieldingPALineshape): - r"""A specialized class for simulating the pure anisotropic spinning sideband - amplitudes of the nuclear shielding resonances resulting from a 2D sideband - separation spectra. - - Args: - anisotropic_dimension: A Dimension object, or an equivalent dictionary - object. This dimension must represent the pure anisotropic - dimension. - inverse_dimension: A list of two Dimension objects, or equivalent - dictionary objects representing the `x`-`y` coordinate grid. - channel: The isotope symbol of the nuclei given as the atomic number - followed by the atomic symbol, for example, `1H`, `13C`, and - `29Si`. This nucleus must correspond to the recorded frequency - resonances. - magnetic_flux_density: The magnetic flux density of the external static - magnetic field. The default value is 9.4 T. - - **Assumption:** - The simulated line-shapes correspond to a finite speed spectrum spinning at the - magic angle, :math:`54.735^\circ`, where the spin rate is the increment along - the anisotropic dimension. - """ - - def __init__( - self, - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density="9.4 T", - ): - - super().__init__( - anisotropic_dimension, - inverse_dimension, - channel, - magnetic_flux_density, - "54.735 deg", - None, - None, - ) - - -# class DAS(LineShape): -# def __init__( -# self, -# anisotropic_dimension, -# inverse_kernel_dimension, -# channel, -# magnetic_flux_density="9.4 T", -# rotor_angle="54.735 deg", -# rotor_frequency="600 Hz", -# number_of_sidebands=None, -# ): -# super().__init__( -# anisotropic_dimension, -# inverse_kernel_dimension, -# channel, -# magnetic_flux_density, -# rotor_angle, -# rotor_frequency, -# number_of_sidebands, -# # "DAS", -# ) - -# def kernel(self, supersampling): -# method = BlochDecayCentralTransitionSpectrum.parse_dict_with_units( -# self.method_args -# ) -# isotope = self.method_args["channels"][0] -# Cq, eta = self._get_zeta_eta(supersampling) -# spin_systems = [ -# SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))]) -# for cq_, e in zip(Cq, eta) -# ] - -# self.simulator.spin_systems = spin_systems -# self.simulator.methods = [method] -# self.simulator.run(pack_as_csdm=False) - -# amp = self.simulator.methods[0].simulation - -# return self._averaged_kernel(amp, supersampling) +from copy import deepcopy + +from mrsimulator import Simulator +from mrsimulator import SpinSystem +from mrsimulator.method.lib import BlochDecaySpectrum + +from mrinversion.kernel.base import LineShape + + +class ShieldingPALineshape(LineShape): + """ + A generalized class for simulating the pure anisotropic NMR nuclear shielding + line-shape kernel. + + Args: + anisotropic_dimension: A Dimension object, or an equivalent dictionary + object. This dimension must represent the pure anisotropic + dimension. + inverse_dimension: A list of two Dimension objects, or equivalent + dictionary objects representing the `x`-`y` coordinate grid. + channel: The channel is an isotope symbol of the nuclei given as the atomic + number followed by the atomic symbol, for example, `1H`, `13C`, and + `29Si`. This nucleus must correspond to the recorded frequency + resonances. + magnetic_flux_density: The magnetic flux density of the external static + magnetic field. The default value is 9.4 T. + rotor_angle: The angle of the sample holder (rotor) relative to the + direction of the external magnetic field. The default value is + 54.735 deg (magic angle). + rotor_frequency: The effective sample spin rate. Depending on the NMR + sequence, this value may be less than the physical sample rotation + frequency. The default is 14 kHz. + number_of_sidebands: The number of sidebands to simulate along the + anisotropic dimension. The default value is 1. + """ + + def __init__( + self, + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="14 kHz", + number_of_sidebands=1, + ): + super().__init__( + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + + def kernel(self, supersampling=1): + """ + Return the NMR nuclear shielding anisotropic line-shape kernel. + + Args: + supersampling: An integer. Each cell is supersampled by the factor + `supersampling` along every dimension. + Returns: + A numpy array containing the line-shape kernel. + """ + args_ = deepcopy(self.method_args) + method = BlochDecaySpectrum.parse_dict_with_units(args_) + isotope = args_["channels"][0] + zeta, eta = self._get_zeta_eta(supersampling) + + x_csdm = self.inverse_kernel_dimension[0] + if x_csdm.coordinates.unit.physical_type == "frequency": + # convert zeta to ppm if given in frequency units. + zeta /= self.larmor_frequency # zeta in ppm + + for dim_i in self.inverse_kernel_dimension: + if dim_i.origin_offset.value == 0: + dim_i.origin_offset = f"{abs(self.larmor_frequency)} MHz" + + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, shielding_symmetric=dict(zeta=z, eta=e))] + ) + for z, e in zip(zeta, eta) + ] + + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + return self._averaged_kernel(amp, supersampling) + + +class MAF(ShieldingPALineshape): + r"""A specialized class for simulating the pure anisotropic NMR nuclear shielding + line-shape kernel resulting from the 2D MAF spectra. + + Args: + anisotropic_dimension: A Dimension object, or an equivalent dictionary + object. This dimension must represent the pure anisotropic + dimension. + inverse_dimension: A list of two Dimension objects, or equivalent + dictionary objects representing the `x`-`y` coordinate grid. + channel: The isotope symbol of the nuclei given as the atomic number + followed by the atomic symbol, for example, `1H`, `13C`, and + `29Si`. This nucleus must correspond to the recorded frequency + resonances. + magnetic_flux_density: The magnetic flux density of the external static + magnetic field. The default value is 9.4 T. + + **Assumptions:** + The simulated line-shapes correspond to an infinite speed spectrum spinning at + :math:`90^\circ`. + """ + + def __init__( + self, + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density="9.4 T", + ): + + super().__init__( + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density, + "90 deg", + "200 GHz", + 1, + ) + + +class SpinningSidebands(ShieldingPALineshape): + r"""A specialized class for simulating the pure anisotropic spinning sideband + amplitudes of the nuclear shielding resonances resulting from a 2D sideband + separation spectra. + + Args: + anisotropic_dimension: A Dimension object, or an equivalent dictionary + object. This dimension must represent the pure anisotropic + dimension. + inverse_dimension: A list of two Dimension objects, or equivalent + dictionary objects representing the `x`-`y` coordinate grid. + channel: The isotope symbol of the nuclei given as the atomic number + followed by the atomic symbol, for example, `1H`, `13C`, and + `29Si`. This nucleus must correspond to the recorded frequency + resonances. + magnetic_flux_density: The magnetic flux density of the external static + magnetic field. The default value is 9.4 T. + + **Assumption:** + The simulated line-shapes correspond to a finite speed spectrum spinning at the + magic angle, :math:`54.735^\circ`, where the spin rate is the increment along + the anisotropic dimension. + """ + + def __init__( + self, + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density="9.4 T", + ): + + super().__init__( + anisotropic_dimension, + inverse_dimension, + channel, + magnetic_flux_density, + "54.735 deg", + None, + None, + ) diff --git a/mrinversion/kernel/nmr.py b/mrinversion/kernel/nmr.py index d3bd1bf..c648566 100644 --- a/mrinversion/kernel/nmr.py +++ b/mrinversion/kernel/nmr.py @@ -1,7 +1,10 @@ -from mrinversion.kernel.csa_aniso import MAF # NOQA -from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA -from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA - -# from mrinversion.kernel.quad_aniso import MQMAS # NOQA - -# from mrinversion.kernel.lineshape import DAS # NOQA +from mrinversion.kernel.csa_aniso import MAF # NOQA +from mrinversion.kernel.csa_aniso import ShieldingPALineshape # NOQA +from mrinversion.kernel.csa_aniso import SpinningSidebands # NOQA +from mrinversion.kernel.quad_aniso import DAS # NOQA +from mrinversion.kernel.quad_aniso import MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMAS # NOQA +from mrinversion.kernel.quad_aniso import SL_MQMASnodist # NOQA + + +# from mrinversion.kernel.quad_aniso import MQMAS # NOQA diff --git a/mrinversion/kernel/quad_aniso.py b/mrinversion/kernel/quad_aniso.py index 2f87bc0..3f85d1f 100644 --- a/mrinversion/kernel/quad_aniso.py +++ b/mrinversion/kernel/quad_aniso.py @@ -1,110 +1,417 @@ -# from copy import deepcopy -# import numpy as np -# from mrsimulator import Simulator -# from mrsimulator import SpinSystem -# from mrsimulator.method import Method -# from mrinversion.kernel.base import LineShape -# class MQMAS(LineShape): -# """ -# A generalized class for simulating the pure anisotropic NMR nuclear shielding -# line-shape kernel. -# Args: -# anisotropic_dimension: A Dimension object, or an equivalent dictionary -# object. This dimension must represent the pure anisotropic -# dimension. -# inverse_dimension: A list of two Dimension objects, or equivalent -# dictionary objects representing the `x`-`y` coordinate grid. -# channel: The channel is an isotope symbol of the nuclei given as the atomic -# number followed by the atomic symbol, for example, `1H`, `13C`, and -# `29Si`. This nucleus must correspond to the recorded frequency -# resonances. -# magnetic_flux_density: The magnetic flux density of the external static -# magnetic field. The default value is 9.4 T. -# """ -# def __init__( -# self, -# anisotropic_dimension, -# inverse_dimension, -# channel, -# magnetic_flux_density="9.4 T", -# ): -# super().__init__( -# anisotropic_dimension, -# inverse_dimension, -# channel, -# magnetic_flux_density, -# rotor_angle=54.7356 * np.pi / 180, -# rotor_frequency=1e9, -# number_of_sidebands=1, -# ) -# def kernel(self, supersampling=1): -# """ -# Return the quadrupolar anisotropic line-shape kernel. -# Args: -# supersampling: An integer. Each cell is supersampled by the factor -# `supersampling` along every dimension. -# Returns: -# A numpy array containing the line-shape kernel. -# """ -# self.inverse_kernel_dimension[0].application["half"] = True -# args_ = deepcopy(self.method_args) -# args_["spectral_dimensions"][0]["events"] = [ -# {"fraction": 27 / 17, "freq_contrib": ["Quad2_0"]}, -# {"fraction": 1, "freq_contrib": ["Quad2_4"]}, -# ] -# method = Method.parse_dict_with_units(args_) -# isotope = args_["channels"][0] -# zeta, eta = self._get_zeta_eta(supersampling) -# # new_size = zeta.size -# # n1, n2 = [item.count for item in self.inverse_kernel_dimension] -# # index = [] -# # for i in range(n2): -# # i_ = i * supersampling -# # for j in range(supersampling): -# # index = np.append(index, np.arange(n1 - i_) + (i_ + j) * n1 + i_) -# # index = np.asarray(index, dtype=int) -# # print(index) -# # zeta = zeta[index] -# # eta = eta[index] -# # larmor frequency from method. -# B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density # in T -# gamma = method.channels[0].gyromagnetic_ratio # in MHz/T -# larmor_frequency = -gamma * B0 # in MHz -# for dim_i in self.inverse_kernel_dimension: -# if dim_i.origin_offset.value == 0: -# dim_i.origin_offset = f"{abs(larmor_frequency)} MHz" -# spin_systems = [ -# SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=z, eta=e))]) -# for z, e in zip(zeta, eta) -# ] -# dim = method.spectral_dimensions[0] -# if dim.origin_offset == 0: -# dim.origin_offset = larmor_frequency * 1e6 # in Hz -# sim = Simulator() -# sim.config.number_of_sidebands = self.number_of_sidebands -# sim.config.decompose_spectrum = "spin_system" -# sim.spin_systems = spin_systems -# sim.methods = [method] -# sim.run(pack_as_csdm=False) -# amp = sim.methods[0].simulation.real -# # amp2 = np.zeros((new_size, amp.shape[1])) -# # amp2 = (amp.T + amp) / 2.0 -# # amp2[index] = amp -# # print(amp2.shape, amp.shape) -# kernel = self._averaged_kernel(amp, supersampling) -# # print(kernel.shape) -# n1, n2 = [item.count for item in self.inverse_kernel_dimension] -# shape = kernel.shape[0] -# kernel.shape = (shape, n1, n2) -# for i in range(n1): -# for j in range(n2): -# if i > j: -# kernel[:, i, j] = 0 -# kernel.shape = (shape, n1 * n2) -# # if not half: -# # return kernel -# index = np.where(kernel.sum(axis=0) != 0)[0] -# self.inverse_kernel_dimension[0].application["index"] = index.tolist() -# return kernel[:, index] -# # -# # return kernel.reshape(shape, n1 * n2) +import warnings + +import csdmpy as cp +import numpy as np +from mrsimulator import Simulator +from mrsimulator import SpinSystem +from mrsimulator.method import Method +from mrsimulator.utils import get_spectral_dimensions + +from mrinversion.kernel.base import LineShape + + +class DAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + # "DAS", + ) + + def kernel(self, supersampling, eta_bound=1): + # update method for DAS spectra events + das_event = dict( + transition_queries=[{"ch1": {"P": [-1], "D": [0]}}], + freq_contrib=["Quad2_4"], + ) + self.method_args["spectral_dimensions"][0]["events"] = [das_event] + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] + ) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " + "xygrid='positive' or xygrid='negative." + ) + + return self._averaged_kernel(amp, supersampling) + + +class MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + + def kernel(self, supersampling, eta_bound=1): + # update method for DAS spectra events + mqmas_events = [ + dict(fraction=-9 / 50, transition_queries=[{"ch1": {"P": [-3], "D": [0]}}]), + dict(fraction=27 / 50, transition_queries=[{"ch1": {"P": [-1], "D": [0]}}]), + ] + + self.method_args["spectral_dimensions"][0]["events"] = mqmas_events + + method = Method.parse_dict_with_units(self.method_args) + isotope = self.method_args["channels"][0] + + Cq, eta, abundances = self._get_zeta_eta( + supersampling, eta_bound, calc_pos=True + ) + + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + sim = Simulator() + sim.config.number_of_sidebands = self.number_of_sidebands + sim.config.decompose_spectrum = "spin_system" + + sim.spin_systems = spin_systems + sim.methods = [method] + sim.run(pack_as_csdm=False) + + amp = sim.methods[0].simulation.real + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " + "xygrid='positive' or xygrid='negative'." + ) + + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMASnodist(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound=1, nQ=3): + import sys + + sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") + # import src.processing as smproc + import src.simulation as smsim + + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + if eta_bound == 1: + Cq, eta = self._get_zeta_eta(supersampling, eta_bound) + else: + Cq, eta, abundances = self._get_zeta_eta(supersampling, eta_bound) + + if eta_bound == 1: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))] + ) + for cq_, e in zip(Cq, eta) + ] + else: + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + + amp = np.asarray( + [ + smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type="c0_c4", + contribs="c0_c4", + return_array=True, + distorted=False, + ) + for mysys in spin_systems + ] + ) + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " + "xygrid='positive' or xygrid='negative." + ) + + return self._averaged_kernel(amp, supersampling) + + +class SL_MQMAS(LineShape): + def __init__( + self, + anisotropic_dimension, + inverse_kernel_dimension, + channel, + exp_dict, + magnetic_flux_density="9.4 T", + rotor_angle="54.735 deg", + rotor_frequency="600 Hz", + number_of_sidebands=None, + ): + super().__init__( + anisotropic_dimension, + inverse_kernel_dimension, + channel, + magnetic_flux_density, + rotor_angle, + rotor_frequency, + number_of_sidebands, + ) + self.exp_dict = exp_dict + self.anisotropic_dimension = anisotropic_dimension + + def kernel(self, supersampling, eta_bound=1, cq_posneg=True, n_quantum=3): + import sys + + sys.path.insert(0, "/home/lexicon2810/github-repos-WSL/mrsmqmas") + # import src.processing as smproc + import src.simulation as smsim + + # import src.fitting as smfit + + isotope = self.method_args["channels"][0] + + Cq, eta, abundances = self._get_zeta_eta( + supersampling, eta_bound, calc_pos=True + ) + + spin_systems = [ + SpinSystem( + sites=[dict(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e))], + abundance=abun, + ) + for cq_, e, abun in zip(Cq, eta, abundances) + ] + + obj = cp.CSDM(dimensions=[self.anisotropic_dimension]) + spec_dim = get_spectral_dimensions(obj) + operators = smsim.build_spin_matrices(nucleus=isotope, n_quantum=n_quantum) + amp = np.asarray( + [ + smsim.simulate_onesite_lineshape( + self.exp_dict, + mysys, + spec_dim[0], + input_type="c0_c4", + contribs="c0_c4", + return_array=True, + distorted=True, + operators=operators, + ) + for mysys in spin_systems + ] + ) + + warnings.warn( + "This kernel is intended to be used with xygrid='mirrored', since we " + "cannot distinguish +Cq from -Cq from an NMR experiment. If you know the" + "sign of Cq from other means, you can restrict the grid using " + "xygrid='positive' or xygrid='negative." + ) + return self._averaged_kernel(amp, supersampling) + + +# from copy import deepcopy +# import numpy as np +# from mrsimulator import Simulator +# from mrsimulator import SpinSystem +# from mrsimulator.method import Method +# from mrinversion.kernel.base import LineShape +# class MQMAS(LineShape): +# """ +# A generalized class for simulating the pure anisotropic NMR nuclear shielding +# line-shape kernel. +# Args: +# anisotropic_dimension: A Dimension object, or an equivalent dictionary +# object. This dimension must represent the pure anisotropic +# dimension. +# inverse_dimension: A list of two Dimension objects, or equivalent +# dictionary objects representing the `x`-`y` coordinate grid. +# channel: The channel is an isotope symbol of the nuclei given as the atomic +# number followed by the atomic symbol, for example, `1H`, `13C`, and +# `29Si`. This nucleus must correspond to the recorded frequency +# resonances. +# magnetic_flux_density: The magnetic flux density of the external static +# magnetic field. The default value is 9.4 T. +# """ +# def __init__( +# self, +# anisotropic_dimension, +# inverse_dimension, +# channel, +# magnetic_flux_density="9.4 T", +# ): +# super().__init__( +# anisotropic_dimension, +# inverse_dimension, +# channel, +# magnetic_flux_density, +# rotor_angle=54.7356 * np.pi / 180, +# rotor_frequency=1e9, +# number_of_sidebands=1, +# ) +# def kernel(self, supersampling=1): +# """ +# Return the quadrupolar anisotropic line-shape kernel. +# Args: +# supersampling: An integer. Each cell is supersampled by the factor +# `supersampling` along every dimension. +# Returns: +# A numpy array containing the line-shape kernel. +# """ +# self.inverse_kernel_dimension[0].application["half"] = True +# args_ = deepcopy(self.method_args) +# args_["spectral_dimensions"][0]["events"] = [ +# {"fraction": 27 / 17, "freq_contrib": ["Quad2_0"]}, +# {"fraction": 1, "freq_contrib": ["Quad2_4"]}, +# ] +# method = Method.parse_dict_with_units(args_) +# isotope = args_["channels"][0] +# zeta, eta = self._get_zeta_eta(supersampling) +# # new_size = zeta.size +# # n1, n2 = [item.count for item in self.inverse_kernel_dimension] +# # index = [] +# # for i in range(n2): +# # i_ = i * supersampling +# # for j in range(supersampling): +# # index = np.append(index, np.arange(n1 - i_) + (i_ + j) * n1 + i_) +# # index = np.asarray(index, dtype=int) +# # print(index) +# # zeta = zeta[index] +# # eta = eta[index] +# # larmor frequency from method. +# B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density # in T +# gamma = method.channels[0].gyromagnetic_ratio # in MHz/T +# larmor_frequency = -gamma * B0 # in MHz +# for dim_i in self.inverse_kernel_dimension: +# if dim_i.origin_offset.value == 0: +# dim_i.origin_offset = f"{abs(larmor_frequency)} MHz" +# spin_systems = [ +# SpinSystem(sites=[dict(isotope=isotope, quadrupolar=dict(Cq=z, eta=e))]) +# for z, e in zip(zeta, eta) +# ] +# dim = method.spectral_dimensions[0] +# if dim.origin_offset == 0: +# dim.origin_offset = larmor_frequency * 1e6 # in Hz +# sim = Simulator() +# sim.config.number_of_sidebands = self.number_of_sidebands +# sim.config.decompose_spectrum = "spin_system" +# sim.spin_systems = spin_systems +# sim.methods = [method] +# sim.run(pack_as_csdm=False) +# amp = sim.methods[0].simulation.real +# # amp2 = np.zeros((new_size, amp.shape[1])) +# # amp2 = (amp.T + amp) / 2.0 +# # amp2[index] = amp +# # print(amp2.shape, amp.shape) +# kernel = self._averaged_kernel(amp, supersampling) +# # print(kernel.shape) +# n1, n2 = [item.count for item in self.inverse_kernel_dimension] +# shape = kernel.shape[0] +# kernel.shape = (shape, n1, n2) +# for i in range(n1): +# for j in range(n2): +# if i > j: +# kernel[:, i, j] = 0 +# kernel.shape = (shape, n1 * n2) +# # if not half: +# # return kernel +# index = np.where(kernel.sum(axis=0) != 0)[0] +# self.inverse_kernel_dimension[0].application["index"] = index.tolist() +# return kernel[:, index] +# # +# # return kernel.reshape(shape, n1 * n2) diff --git a/mrinversion/kernel/utils.py b/mrinversion/kernel/utils.py index d60ee19..fea9ace 100644 --- a/mrinversion/kernel/utils.py +++ b/mrinversion/kernel/utils.py @@ -1,173 +1,262 @@ -import numpy as np - - -def x_y_to_zeta_eta(x, y): - r"""Convert the coordinates :math:`(x,y)` to :math:`(\zeta, \eta)` using the - following definition, - - .. math:: - \left.\begin{array}{rl} - \zeta &= \sqrt{x^2 + y^2}, \\ - \eta &= \frac{4}{\pi} \tan^{-1} \left| \frac{x}{y} \right| - \end{array} {~~~~~~~~} \right\} {~~~~~~~~} |x| \le |y|. - - .. math:: - \left.\begin{array}{rl} - \zeta &= -\sqrt{x^2 + y^2}, \\ - \eta &= \frac{4}{\pi} \tan^{-1} \left| \frac{y}{x} \right| - \end{array} {~~~~~~~~} \right\} {~~~~~~~~} |x| > |y|. - - Args: - x: floats or Quantity object. The coordinate x. - y: floats or Quantity object. The coordinate y. - - Return: - A list of two ndarrays. The first array is the :math:`\zeta` - coordinates. The second array is the :math:`\eta` coordinates. - """ - x_unit = y_unit = 1 - if x.__class__.__name__ == "Quantity": - x_unit = x.unit - x = x.value - if y.__class__.__name__ == "Quantity": - y_unit = y.unit - y = y.value - if x_unit != y_unit: - raise ValueError( - f"x and y must have same dimensionality; x ({x_unit}) != y ({y_unit})." - ) - - zeta = np.sqrt(x**2 + y**2) # + offset - eta = 1.0 - if x > y: - zeta = -zeta - eta = (4.0 / np.pi) * np.arctan(y / x) - - if x < y: - eta = (4.0 / np.pi) * np.arctan(x / y) - - return zeta * x_unit, eta - - -def _x_y_to_zeta_eta(x, y): - """Same as def x_y_to_zeta_eta, but for ndarrays.""" - x = np.abs(x) - y = np.abs(y) - zeta = np.sqrt(x**2 + y**2) # + offset - eta = np.ones(zeta.shape) - index = np.where(x > y) - zeta[index] = -zeta[index] - eta[index] = (4.0 / np.pi) * np.arctan(y[index] / x[index]) - - index = np.where(x < y) - eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index]) - - return zeta.ravel(), eta.ravel() - - -def zeta_eta_to_x_y(zeta, eta): - r"""Convert the coordinates :math:`(\zeta,\eta)` to :math:`(x, y)` using the - following definition, - - .. math:: - \left. \begin{array}{rl} - x &= |\zeta| \sin\theta, \\ - y &= |\zeta| \cos\theta - \end{array} {~~~~~~~~} \right\} {~~~~~~~~} \zeta \ge 0 - - .. math:: - \left. \begin{array}{rl} - x &= |\zeta| \cos\theta, \\ - y &= |\zeta| \sin\theta - \end{array} {~~~~~~~~} \right\} {~~~~~~~~} \zeta < 0 - - where :math:`\theta = \frac{\pi}{4}\eta`. - - Args: - x: ndarray or list of floats. The coordinate x. - y: ndarray or list of floats. The coordinate y. - - Return: - A list of ndarrays. The first array holds the coordinate :math:`x`. The - second array holds the coordinates :math:`y`. - """ - zeta = np.asarray(zeta) - eta = np.asarray(eta) - - theta = np.pi * eta / 4.0 - x = np.zeros(zeta.size) - y = np.zeros(zeta.size) - - index = np.where(zeta >= 0) - x[index] = zeta[index] * np.sin(theta[index]) - y[index] = zeta[index] * np.cos(theta[index]) - - index = np.where(zeta < 0) - x[index] = -zeta[index] * np.cos(theta[index]) - y[index] = -zeta[index] * np.sin(theta[index]) - - return x.ravel(), y.ravel() - - -def _x_y_to_zeta_eta_distribution(grid, supersampling): - """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" - x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) - y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) - - if x_coordinates.unit.physical_type == "frequency": - x_coordinates = x_coordinates.to("Hz").value - y_coordinates = y_coordinates.to("Hz").value - - elif x_coordinates.unit.physical_type == "dimensionless": - x_coordinates = x_coordinates.to("ppm").value - y_coordinates = y_coordinates.to("ppm").value - - x_mesh, y_mesh = np.meshgrid( - np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" - ) - - return _x_y_to_zeta_eta(x_mesh, y_mesh) - - -def _supersampled_coordinates(dimension, supersampling=1): - r"""The coordinates along the dimension. - - Args: - supersampling: An integer used in supersampling the coordinates along the - dimension, If :math:`n` is the count, :math:`\Delta_x` is the increment, - :math:`x_0` is the coordinates offset along the dimension, and :math:`m` is - the supersampling factor, a total of :math:`mn` coordinates are sampled - using - - .. math:: - x = [0 .. (nm-1)] \Delta_x + \x_0 - \frac{1}{2} \Delta_x (m-1) - - where :math:`\Delta_x' = \frac{\Delta_x}{m}`. - - Returns: - An `Quantity` array of coordinates. - """ - if dimension.type == "linear": - increment = dimension.increment / supersampling - array = np.arange(dimension.count * supersampling) * increment - array += dimension.coordinates_offset - # shift the coordinates by half a bin for proper averaging - array -= 0.5 * increment * (supersampling - 1) - - if dimension.type == "monotonic": - coordinates = dimension.coordinates - unit = coordinates[0].unit - size = coordinates.size - - diff = np.zeros(size) - diff[1:] = (coordinates[1:] - coordinates[:-1]) / supersampling - diff *= unit - - s2 = supersampling // 2 - eo = 0.5 if supersampling % 2 == 0 else 0 - - array = np.zeros(size * supersampling) * unit - for i in range(supersampling): - array[i::supersampling] = coordinates + (i - s2 + eo) * diff - - return array +import numpy as np + + +def x_y_to_zeta_eta(x, y): + r"""Convert the coordinates :math:`(x,y)` to :math:`(\zeta, \eta)` using the + following definition, + + .. math:: + \left.\begin{array}{rl} + \zeta &= \sqrt{x^2 + y^2}, \\ + \eta &= \frac{4}{\pi} \tan^{-1} \left| \frac{x}{y} \right| + \end{array} {~~~~~~~~} \right\} {~~~~~~~~} |x| \le |y|. + + .. math:: + \left.\begin{array}{rl} + \zeta &= -\sqrt{x^2 + y^2}, \\ + \eta &= \frac{4}{\pi} \tan^{-1} \left| \frac{y}{x} \right| + \end{array} {~~~~~~~~} \right\} {~~~~~~~~} |x| > |y|. + + Args: + x: floats or Quantity object. The coordinate x. + y: floats or Quantity object. The coordinate y. + + Return: + A list of two ndarrays. The first array is the :math:`\zeta` + coordinates. The second array is the :math:`\eta` coordinates. + """ + x_unit = y_unit = 1 + if x.__class__.__name__ == "Quantity": + x_unit = x.unit + x = x.value + if y.__class__.__name__ == "Quantity": + y_unit = y.unit + y = y.value + if x_unit != y_unit: + raise ValueError( + f"x and y must have same dimensionality: x ({x_unit}) != y ({y_unit})." + ) + + zeta = np.sqrt(x**2 + y**2) # + offset + eta = 1.0 + if x > y: + zeta = -zeta + eta = (4.0 / np.pi) * np.arctan(y / x) + + if x < y: + eta = (4.0 / np.pi) * np.arctan(x / y) + + return zeta * x_unit, eta + + +def _x_y_to_zeta_eta(x, y, eta_bound=1, calc_pos=False): + """Same as def x_y_to_zeta_eta, but for ndarrays.""" + x = np.abs(x).ravel() + y = np.abs(y).ravel() + zeta = np.sqrt(x**2 + y**2) # + offset + abundances = np.ones(zeta.shape) + eta = np.ones(zeta.shape) + eta = eta.tolist() + zeta = zeta.tolist() + del_these = [] + for index, _ in enumerate(x): + if x[index] > y[index]: + if not calc_pos: + zeta[index] = -zeta[index] + this_eta = (4.0 / np.pi) * np.arctan(y[index] / x[index]) + if this_eta < eta_bound: + eta[index] = this_eta + else: + abundances[index] = 0 + else: + abundances[index] = 0 + + elif x[index] < y[index]: + this_eta = (4.0 / np.pi) * np.arctan(x[index] / y[index]) + if this_eta < eta_bound: + eta[index] = this_eta + else: + abundances[index] = 0 + elif x[index] == y[index] and eta_bound < 1: + abundances[index] = 0 + + if eta_bound == 1 and not calc_pos: + return np.asarray(zeta), np.asarray(eta) + else: + for idx in del_these[::-1]: + del zeta[idx], eta[idx] + return np.asarray(zeta), np.asarray(eta), abundances + + +def zeta_eta_to_x_y(zeta, eta): + r"""Convert the coordinates :math:`(\zeta,\eta)` to :math:`(x, y)` using the + following definition, + + .. math:: + \left. \begin{array}{rl} + x &= |\zeta| \sin\theta, \\ + y &= |\zeta| \cos\theta + \end{array} {~~~~~~~~} \right\} {~~~~~~~~} \zeta \ge 0 + + .. math:: + \left. \begin{array}{rl} + x &= |\zeta| \cos\theta, \\ + y &= |\zeta| \sin\theta + \end{array} {~~~~~~~~} \right\} {~~~~~~~~} \zeta < 0 + + where :math:`\theta = \frac{\pi}{4}\eta`. + + Args: + x: ndarray or list of floats. The coordinate x. + y: ndarray or list of floats. The coordinate y. + + Return: + A list of ndarrays. The first array holds the coordinate :math:`x`. The + second array holds the coordinates :math:`y`. + """ + zeta = np.asarray(zeta) + eta = np.asarray(eta) + + theta = np.pi * eta / 4.0 + x = np.zeros(zeta.size) + y = np.zeros(zeta.size) + + index = np.where(zeta >= 0) + x[index] = zeta[index] * np.sin(theta[index]) + y[index] = zeta[index] * np.cos(theta[index]) + + index = np.where(zeta < 0) + x[index] = -zeta[index] * np.cos(theta[index]) + y[index] = -zeta[index] * np.sin(theta[index]) + + return x.ravel(), y.ravel() + + +def _x_y_to_zeta_eta_distribution(grid, supersampling, eta_bound=1, calc_pos=False): + """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" + x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) + y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) + + if x_coordinates.unit.physical_type == "frequency": + x_coordinates = x_coordinates.to("Hz").value + y_coordinates = y_coordinates.to("Hz").value + + elif x_coordinates.unit.physical_type == "dimensionless": + x_coordinates = x_coordinates.to("ppm").value + y_coordinates = y_coordinates.to("ppm").value + + x_mesh, y_mesh = np.meshgrid( + np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" + ) + + return _x_y_to_zeta_eta(x_mesh, y_mesh, eta_bound, calc_pos) + + +# def _x_y_to_cq_eta_distribution(grid, supersampling): +# """Return a list of zeta-eta coordinates from a list of x-y coordinates.""" +# x_coordinates = _supersampled_coordinates(grid[0], supersampling=supersampling) +# y_coordinates = _supersampled_coordinates(grid[1], supersampling=supersampling) + +# if x_coordinates.unit.physical_type == "frequency": +# x_coordinates = x_coordinates.to("Hz").value +# y_coordinates = y_coordinates.to("Hz").value + +# elif x_coordinates.unit.physical_type == "dimensionless": +# x_coordinates = x_coordinates.to("ppm").value +# y_coordinates = y_coordinates.to("ppm").value + +# x_mesh, y_mesh = np.meshgrid( +# np.abs(x_coordinates), np.abs(y_coordinates), indexing="xy" +# ) + +# return _x_y_to_cq_eta(x_mesh, y_mesh) + + +def _supersampled_coordinates(dimension, supersampling=1): + r"""The coordinates along the dimension. + + Args: + supersampling: An integer used in supersampling the coordinates along the + dimension, If :math:`n` is the count, :math:`\Delta_x` is the increment, + :math:`x_0` is the coordinates offset along the dimension, and :math:`m` is + the supersampling factor, a total of :math:`mn` coordinates are sampled + using + + .. math:: + x = [0 .. (nm-1)] \Delta_x + \x_0 - \frac{1}{2} \Delta_x (m-1) + + where :math:`\Delta_x' = \frac{\Delta_x}{m}`. + + Returns: + An `Quantity` array of coordinates. + """ + if dimension.type == "linear": + increment = dimension.increment / supersampling + array = np.arange(dimension.count * supersampling) * increment + array += dimension.coordinates_offset + # shift the coordinates by half a bin for proper averaging + array -= 0.5 * increment * (supersampling - 1) + + if dimension.type == "monotonic": + coordinates = dimension.coordinates + unit = coordinates[0].unit + size = coordinates.size + + diff = np.zeros(size) + diff[1:] = (coordinates[1:] - coordinates[:-1]) / supersampling + diff *= unit + + s2 = supersampling // 2 + eo = 0.5 if supersampling % 2 == 0 else 0 + + array = np.zeros(size * supersampling) * unit + for i in range(supersampling): + array[i::supersampling] = coordinates + (i - s2 + eo) * diff + + return array + + +# def cq_eta_to_x_y(cq, eta): +# r"""Convert the coordinates :math:`(C_q,\eta)` to :math:`(x, y)` using the +# following definition, + +# .. math:: +# \begin{array}{rl} +# x &= C_q \sin\theta, \\ +# y &= C_q \cos\theta +# \end{array} {~~~~~~~~} + +# where :math:`\theta = \frac{\pi}{2}\eta`. + +# Args: +# x: ndarray or list of floats. The coordinate x. +# y: ndarray or list of floats. The coordinate y. + +# Return: +# A list of ndarrays. The first array holds the coordinate :math:`x`. The +# second array holds the coordinates :math:`y`. +# """ +# cq = np.asarray(cq) +# eta = np.asarray(eta) + +# theta = np.pi * eta / 2.0 +# x = np.zeros(cq.size) +# y = np.zeros(cq.size) + +# index = np.arange(len(cq)) +# x[index] = cq[index] * np.cos(theta[index]) +# y[index] = cq[index] * np.sin(theta[index]) + +# return x.ravel(), y.ravel() + + +# def _x_y_to_cq_eta(x, y): +# """Same as def x_y_to_zeta_eta, but for ndarrays.""" +# x = np.abs(x) +# y = np.abs(y) +# cq = np.sqrt(x**2 + y**2) # + offset +# eta = np.ones(cq.shape) +# index = np.arange(len(cq)) +# eta[index] = (2.0 / np.pi) * np.arctan(y[index] / x[index]) + +# return cq.ravel(), eta.ravel() diff --git a/mrinversion/linear_model/_base_l1l2.py b/mrinversion/linear_model/_base_l1l2.py index 27b9ed1..65c8932 100644 --- a/mrinversion/linear_model/_base_l1l2.py +++ b/mrinversion/linear_model/_base_l1l2.py @@ -1,677 +1,765 @@ -from copy import deepcopy - -import csdmpy as cp -import matplotlib.pyplot as plt -import numpy as np -from joblib import delayed -from joblib import Parallel -from sklearn.linear_model import Lasso -from sklearn.linear_model import LassoLars -from sklearn.linear_model import MultiTaskLasso -from sklearn.model_selection import cross_validate -from sklearn.model_selection import KFold - - -__author__ = "Deepansh J. Srivastava" -__email__ = "srivastava.89@osu.edu" - - -class GeneralL2Lasso: - r"""The Minimizer class solves the following equation, - - .. math:: - {\bf f} = \underset{{\bf f}}{\text{argmin}} \left( \frac{1}{m} \| - {\bf Kf - s} \|^2_2 + - \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 + - \lambda \| {\bf f} \|_1 \right), - - where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, - :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known signal - containing noise, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` - is the desired solution matrix. - - - Based on the regularization literal, the above problem is constraint - - Args: - alpha: Float, the hyperparameter, :math:`\alpha`. - lambda1: Float, the hyperparameter, :math:`\lambda`. - max_iterations: Integer, the maximum number of iterations allowed when - solving the problem. The default value is 10000. - tolerance: Float, the tolerance at which the solution is - considered converged. The default value is 1e-5. - positive: Boolean. If True, the amplitudes in the solution, - :math:`{\bf f}` is all positive, else the solution may contain - positive and negative amplitudes. The default is True. - regularizer: String, a literal specifying the form of matrix - :math:`{\bf J}_i`. The allowed literals are `smooth lasso` - and `sparse ridge fusion`. - f_shape: The shape of the solution, :math:`{\bf f}`, given as a tuple - (n1, n2, ..., nd) - """ - - 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 - - def fit(self, K, s): - r"""Fit the model using the coordinate descent method from scikit-learn. - - Args - ---- - - K: ndarray - The :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of - shape (m, n). - s: ndarray or CSDM object. - A csdm object or an equivalent numpy array holding the signal, - :math:`{\bf s}`, as a :math:`m \times m_\text{count}` matrix. - """ - s_, self.scale = prepare_signal(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}." - ) - - alpha = s_.size * self.hyperparameters["alpha"] - Ks, ss = _get_augmented_data( - K=K, s=s_, alpha=alpha, regularizer=self.regularizer, f_shape=self.f_shape - ) - - # 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": - estimator = MultiTaskLasso( - alpha=self.hyperparameters["lambda"] / 2.0, - fit_intercept=False, - copy_X=True, - max_iter=self.max_iterations, - tol=self.tolerance, - warm_start=False, - random_state=None, - selection="random", - # positive=self.positive, - ) - - if self.method == "gradient_decent": - estimator = Lasso( - alpha=self.hyperparameters["lambda"] / 2.0, - fit_intercept=False, - copy_X=True, - max_iter=self.max_iterations, - tol=self.tolerance, - warm_start=False, - random_state=None, - selection="random", - positive=self.positive, - ) - - if self.method == "lars": - estimator = LassoLars( - alpha=self.hyperparameters["lambda"] / 2.0, - fit_intercept=False, - verbose=True, - # normalize=False, - precompute=True, - max_iter=self.max_iterations, - eps=2.220446049250313e-16, - copy_X=True, - fit_path=False, - positive=True, - jitter=None, - random_state=None, - ) - - estimator.fit(Ks, ss) - f = estimator.coef_.copy() - if s_.shape[1] > 1 and len(self.f_shape) == 2: - f.shape = (s_.shape[1],) + self.f_shape - f[:, :, 0] /= 2.0 - f[:, 0, :] /= 2.0 - elif s_.shape[1] == 1 and len(self.f_shape) == 2: - f.shape = self.f_shape - f[:, 0] /= 2.0 - f[0, :] /= 2.0 - - f *= self.scale - self.estimator = estimator - self.f = f - self.n_iter = estimator.n_iter_ - self._sol_to_csdm(s) - - def _sol_to_csdm(self, s): - """Pack solution as csdm object""" - if not isinstance(s, cp.CSDM): - return - - self.f = cp.as_csdm(self.f) - - for i, dim in enumerate(self.inverse_dimension): - app = dim.application or {} - if "com.github.deepanshs.mrinversion" in app: - meta = app["com.github.deepanshs.mrinversion"] - is_log = meta.get("log", False) - if is_log: - coords = np.log10(dim.coordinates.value) - self.inverse_dimension[i] = cp.as_dimension( - array=coords, label=meta["label"] - ) - - if len(s.dimensions) > 1 and len(self.f.shape) == 3: - self.f.dimensions[2] = s.dimensions[1] - self.f.dimensions[1] = self.inverse_dimension[1] - self.f.dimensions[0] = self.inverse_dimension[0] - - elif len(s.dimensions) == 1 and len(self.f.shape) == 2: - self.f.dimensions[1] = self.inverse_dimension[1] - self.f.dimensions[0] = self.inverse_dimension[0] - - elif len(s.dimensions) > 1: - self.f.dimensions[1] = s.dimensions[1] - self.f.dimensions[0] = self.inverse_dimension[0] - - 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 - """ - predict = self.estimator.predict(K) * self.scale - return predict - - def residuals(self, K, s): - r"""Return the residual as the difference the data and the predicted 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.copy()) - 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 GeneralL2LassoCV: - def __init__( - self, - alphas, - lambdas, - 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.cv_alphas = np.asarray(alphas).ravel() - 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] - - def fit(self, K, s): - 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. - - 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). - """ - s_, self.scale = prepare_signal(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}." - ) - - 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]) - - Ks, ss = _get_augmented_data( - K=K, - s=s_, - alpha=s_.size * self.cv_alphas[0], - regularizer=self.regularizer, - f_shape=self.f_shape, - ) - start_index = K.shape[0] - - l1 = self._get_minimizer() - # l1.fit(Ks, ss) - - l1_array = [] - for lambda_ in self.cv_lambdas: - l1_array.append(deepcopy(l1)) - l1_array[-1].alpha = lambda_ / 2.0 - - j = 0 - for alpha_ratio_ in alpha_ratio: - if alpha_ratio_ != 0: - Ks[start_index:] *= alpha_ratio_ - jobs = ( - delayed(cv)(l1_array[i], Ks, ss, cv_indexes) - for i in range(self.cv_lambdas.size) - ) - 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. - self.cv_map -= (self.sigma / self.scale) ** 2 - - # 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. - 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( - alpha=self.hyperparameters["alpha"], - lambda1=self.hyperparameters["lambda"], - max_iterations=self.max_iterations, - tolerance=self.tolerance, - positive=self.positive, - regularizer=self.regularizer, - inverse_dimension=self.inverse_dimension, - method=self.method, - ) - self.opt.fit(K, s) - self.f = self.opt.f - self.cv_map_to_csdm() - - def cv_map_to_csdm(self): - # convert cv_map to csdm - self.cv_map = cp.as_csdm(np.squeeze(self.cv_map.T.copy())) - 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: - return - - d1 = cp.as_dimension(-np.log10(self.cv_lambdas), label="-log(Îģ)") - if len(self.cv_alphas) != 1: - self.cv_map.dimensions[1] = d1 - else: - self.cv_map.dimensions[0] = d1 - - def cv_plot(self): - plt.figure(figsize=(5, 3.5)) - ax = plt.subplot(projection="csdm") - ax.contour(np.log10(self.cv_map), levels=25) - ax.scatter( - -np.log10(self.hyperparameters["alpha"]), - -np.log10(self.hyperparameters["lambda"]), - marker="x", - color="k", - ) - plt.tight_layout(pad=0.5) - plt.show() - - def _get_minimizer(self): - """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=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=True, - selection="random", - ) - - if self.method == "gradient_decent": - return Lasso( - alpha=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=True, - selection="random", - ) - - if self.method == "lars": - return LassoLars( - alpha=self.cv_lambdas[0] / 2.0, - fit_intercept=False, - verbose=True, - # normalize=False, - precompute="auto", - max_iter=self.max_iterations, - eps=2.220446049250313e-16, - copy_X=True, - fit_path=False, - positive=self.positive, - jitter=None, - random_state=None, - ) - - 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). - - 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 predicted 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) - - @property - def cross_validation_curve(self): - """The cross-validation error metric determined as the mean square error. - - Returns: A two-dimensional CSDM object.""" - return self.cv_map - - -def cv(l1, X, y, cv): - """Return the cross-validation score as negative of mean square error.""" - # if isinstance(l1, (Lasso, MultiTaskLasso)): - # fit_params = {"check_input": False} - # if isinstance(l1, LassoLars): - # fit_params = None # {"Xy": np.dot(X.T, y)} - - cv_score = cross_validate( - l1, - X=X, - y=y, - scoring="neg_mean_squared_error", # 'neg_mean_absolute_error", - cv=cv, - # fit_params=fit_params, - n_jobs=1, - verbose=0, - ) - return cv_score["test_score"].mean() - - -def _get_smooth_size(f_shape, regularizer, max_size): - r"""Return the number of rows appended to for the augmented kernel. - - For smooth-lasso, the number of rows is given as - rows = \prod_{i=1}^d n_i (\sum_{j=0}^d (n_j-1)/n_j) - - For sparse ridge fusion, the number of rows is given as - rows = \prod_{i=1}^d n_i (\sum_{j=0}^d (n_j-2)/n_j) - """ - shape = np.asarray(f_shape) - shape_prod = shape.prod() - if shape_prod != max_size: - raise ValueError( - "The product of the shape must be equal to the length of axis 1 of the " - "kernel, K" - ) - if regularizer == "smooth lasso": - smooth_array_size = [int(shape_prod * (i - 1) / i) for i in shape] - smooth_size = np.asarray(smooth_array_size).sum() - elif regularizer == "sparse ridge fusion": - smooth_array_size = [int(shape_prod * (i - 2) / i) for i in shape] - smooth_size = np.asarray(smooth_array_size).sum() - else: - smooth_size = 0 - - return smooth_size - - -def _get_cv_indexes(K, folds, regularizer, f_shape=None, random=False, times=1): - """Return the indexes of the kernel and signal, corresponding to the test - and train sets.""" - cv_indexes = [] - ks0, ks1 = K.shape - if isinstance(f_shape, int): - f_shape = (f_shape,) - - tr_ = [] - smooth_size = _get_smooth_size(f_shape, regularizer, ks1) - - tr_ = (np.arange(smooth_size) + ks0).tolist() - - for j in range(folds): - train_ = [] - test_ = [] - for i in range(ks0): - if (i + j - folds) % folds == 0: - test_.append(i) - else: - train_.append(i) - train_ += tr_ - cv_indexes.append([train_, test_]) - - if random: - for _ in range(times): - kf = KFold(n_splits=folds, shuffle=True) - kf.get_n_splits(K) - for train_index, test_index in kf.split(K): - cv_indexes.append([train_index.tolist() + tr_, test_index]) - return cv_indexes - - -def generate_J_i(Ai, alpha, f_shape): - J = [] - sqrt_alpha = np.sqrt(alpha) - identity = [np.eye(i) for i in f_shape] - for i, i_count in enumerate(f_shape): - J_array = deepcopy(identity) - J_array[i] = Ai(i_count) - Ji_ = 1 - for j_array_ in J_array: - Ji_ = np.kron(Ji_, j_array_) - J.append(Ji_ * sqrt_alpha) - return J - - -def _get_augmented_data(K, s, alpha, regularizer, f_shape=None): - """Creates a smooth kernel, K, with alpha regularization parameter.""" - if alpha == 0: - return np.asfortranarray(K), np.asfortranarray(s) - - ks0, ks1 = K.shape - ss0, ss1 = s.shape - - if isinstance(f_shape, int): - f_shape = (f_shape,) - - smooth_size = _get_smooth_size(f_shape, regularizer, ks1) - - if regularizer == "smooth lasso": - - def Ai_smooth_lasso(i): - return (-1 * np.eye(i) + np.diag(np.ones(i - 1), k=-1))[1:] - - J = generate_J_i(Ai_smooth_lasso, alpha, f_shape) - - if regularizer == "sparse ridge fusion": - - def Ai_sparse_ridge_fusion(i): - A_temp = -1 * np.eye(i) - A_temp += 2 * np.diag(np.ones(i - 1), k=-1) - A_temp += -1 * np.diag(np.ones(i - 2), k=-2) - return A_temp[2:] - # return ( - # -1 * np.eye(i) - # + 2 * np.diag(np.ones(i - 1), k=-1) - # - 1 * np.diag(np.ones(i - 2), k=-2) - # )[2:] - - J = generate_J_i(Ai_sparse_ridge_fusion, alpha, f_shape) - - K_ = np.empty((ks0 + smooth_size, ks1)) - - K_[:ks0] = K - start = ks0 - end = ks0 - for J_i in J: - end = end + J_i.shape[0] - K_[start:end] = J_i - start = end - - s_ = np.zeros((ss0 + smooth_size, ss1)) - s_[:ss0] = s.real - - return np.asfortranarray(K_), np.asfortranarray(s_) - - -def prepare_signal(sig): - """Prepare signal for inversion""" - s_ = sig.y[0].components[0].T if isinstance(sig, cp.CSDM) else sig - s_ = s_[:, np.newaxis] if s_.ndim == 1 else s_ - - scale = np.sqrt(np.mean(np.abs(s_) ** 2)) - s_ = s_ / scale - return s_, scale +from copy import deepcopy + +import csdmpy as cp +import matplotlib.pyplot as plt +import numpy as np +from joblib import delayed +from joblib import Parallel +from sklearn.linear_model import Lasso +from sklearn.linear_model import LassoLars +from sklearn.linear_model import MultiTaskLasso +from sklearn.model_selection import cross_validate +from sklearn.model_selection import KFold + + +__author__ = "Deepansh J. Srivastava" +__email__ = "srivastava.89@osu.edu" +ALLOWED_GRIDS = ("full", "mirrored", "positive", "negative") + + +class GeneralL2Lasso: + r"""The Minimizer class solves the following equation, + + .. math:: + {\bf f} = \underset{{\bf f}}{\text{argmin}} \left( \frac{1}{m} \| + {\bf Kf - s} \|^2_2 + + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 + + \lambda \| {\bf f} \|_1 \right), + + where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, + :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known signal + containing noise, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` + is the desired solution matrix. + + + Based on the regularization literal, the above problem is constraint + + Args: + alpha: Float, the hyperparameter, :math:`\alpha`. + lambda1: Float, the hyperparameter, :math:`\lambda`. + max_iterations: Integer, the maximum number of iterations allowed when + solving the problem. The default value is 10000. + tolerance: Float, the tolerance at which the solution is + considered converged. The default value is 1e-5. + positive: Boolean. If True, the amplitudes in the solution, + :math:`{\bf f}` is all positive, else the solution may contain + positive and negative amplitudes. The default is True. + regularizer: String, a literal specifying the form of matrix + :math:`{\bf J}_i`. The allowed literals are `smooth lasso` + and `sparse ridge fusion`. + f_shape: The shape of the solution, :math:`{\bf f}`, given as a tuple + (n1, n2, ..., nd) + """ + + 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", + xygrid="full", + ): + if xygrid not in ALLOWED_GRIDS: + raise KeyError( + "Please choose a valid option for 'xygrid'. Valid options are 'full', " + "'mirrored', 'positive', and 'negative'." + ) + + self.hyperparameters = {"lambda": lambda1, "alpha": alpha} + self.max_iterations = max_iterations + self.tolerance = tolerance + 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 + self.xygrid = xygrid + + def fit(self, K, s): + r"""Fit the model using the coordinate descent method from scikit-learn. + + Args + ---- + + K: ndarray + The :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of + shape (m, n). + s: ndarray or CSDM object. + A csdm object or an equivalent numpy array holding the signal, + :math:`{\bf s}`, as a :math:`m \times m_\text{count}` matrix. + """ + s_, self.scale = prepare_signal(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}." + ) + + alpha = s_.size * self.hyperparameters["alpha"] + Ks, ss = _get_augmented_data( + K=K, s=s_, alpha=alpha, regularizer=self.regularizer, f_shape=self.f_shape + ) + + # 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": + estimator = MultiTaskLasso( + alpha=self.hyperparameters["lambda"] / 2.0, + fit_intercept=False, + copy_X=True, + max_iter=self.max_iterations, + tol=self.tolerance, + warm_start=False, + random_state=None, + selection="random", + # positive=self.positive, + ) + + if self.method == "gradient_decent": + estimator = Lasso( + alpha=self.hyperparameters["lambda"] / 2.0, + fit_intercept=False, + copy_X=True, + max_iter=self.max_iterations, + tol=self.tolerance, + warm_start=False, + random_state=None, + selection="random", + positive=self.positive, + ) + + if self.method == "lars": + estimator = LassoLars( + alpha=self.hyperparameters["lambda"] / 2.0, + fit_intercept=False, + verbose=True, + # normalize=False, + precompute=True, + max_iter=self.max_iterations, + eps=2.220446049250313e-16, + copy_X=True, + fit_path=False, + positive=True, + jitter=None, + random_state=None, + ) + + estimator.fit(Ks, ss) + f = estimator.coef_.copy() + if s_.shape[1] > 1 and len(self.f_shape) == 2: + f.shape = (s_.shape[1],) + self.f_shape + f[:, :, 0] /= 2.0 + f[:, 0, :] /= 2.0 + elif s_.shape[1] == 1 and len(self.f_shape) == 2: + f.shape = self.f_shape + f[:, 0] /= 2.0 + f[0, :] /= 2.0 + + f *= self.scale + self.estimator = estimator + self.f = f + self.n_iter = estimator.n_iter_ + self._sol_to_csdm(s) + + def _sol_to_csdm(self, s): + """Pack solution as csdm object""" + if not isinstance(s, cp.CSDM): + return + + self.f = cp.as_csdm(self.f) + + for i, dim in enumerate(self.inverse_dimension): + app = dim.application or {} + if "com.github.deepanshs.mrinversion" in app: + meta = app["com.github.deepanshs.mrinversion"] + is_log = meta.get("log", False) + if is_log: + coords = np.log10(dim.coordinates.value) + self.inverse_dimension[i] = cp.as_dimension( + array=coords, label=meta["label"] + ) + + if len(s.dimensions) > 1 and len(self.f.shape) == 3: + self.f.dimensions[2] = s.dimensions[1] + self.f.dimensions[1] = self.inverse_dimension[1] + self.f.dimensions[0] = self.inverse_dimension[0] + + elif len(s.dimensions) == 1 and len(self.f.shape) == 2: + self.f.dimensions[1] = self.inverse_dimension[1] + self.f.dimensions[0] = self.inverse_dimension[0] + + elif len(s.dimensions) > 1: + self.f.dimensions[1] = s.dimensions[1] + self.f.dimensions[0] = self.inverse_dimension[0] + + if self.xygrid == "mirrored": + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() + + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + y[dim0, dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def predict(self, K): + r"""Predict the signal using the linear model. + + 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 + """ + predict = self.estimator.predict(K) * self.scale + return predict + + def residuals(self, K, s): + r"""Return the residual as the difference the data and the predicted 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.copy()) + 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 GeneralL2LassoCV: + def __init__( + self, + alphas, + lambdas, + 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", + xygrid="full", + ): + if xygrid not in ALLOWED_GRIDS: + raise KeyError( + "Please choose a valid option for 'xygrid'. Valid options are 'full', " + "'mirrored', 'positive', and 'negative'." + ) + + self.cv_alphas = np.asarray(alphas).ravel() + self.cv_lambdas = np.asarray(lambdas).ravel() + + 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.xygrid = xygrid + + def fit(self, K, s, cv_map_as_csdm=True): + r"""Fit the model using the coordinate descent method from scikit-learn for + all alpha anf lambda values using the `n`-folds cross-validation technique. + The cross-validation metric is the mean squared error. + + 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). + """ + s_, self.scale = prepare_signal(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}." + ) + + 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]) + + Ks, ss = _get_augmented_data( + K=K, + s=s_, + alpha=s_.size * self.cv_alphas[0], + regularizer=self.regularizer, + f_shape=self.f_shape, + ) + start_index = K.shape[0] + + l1 = self._get_minimizer() + # l1.fit(Ks, ss) + + l1_array = [] + for lambda_ in self.cv_lambdas: + l1_array.append(deepcopy(l1)) + l1_array[-1].alpha = lambda_ / 2.0 + + j = 0 + for alpha_ratio_ in alpha_ratio: + if alpha_ratio_ != 0: + Ks[start_index:] *= alpha_ratio_ + jobs = ( + delayed(cv)(l1_array[i], Ks, ss, cv_indexes) + for i in range(self.cv_lambdas.size) + ) + 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. + self.cv_map -= (self.sigma / self.scale) ** 2 + + # 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. + 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( + alpha=self.hyperparameters["alpha"], + lambda1=self.hyperparameters["lambda"], + max_iterations=self.max_iterations, + tolerance=self.tolerance, + positive=self.positive, + regularizer=self.regularizer, + inverse_dimension=self.inverse_dimension, + method=self.method, + xygrid=self.xygrid, + ) + + self.opt.fit(K, s) + self.f = self.opt.f + + if self.xygrid == "mirrored": + # print('made it to mirrored') + self.f = self.mirror_inversion() + elif self.xygrid == "negative": + self.f = self.flip_inversion() + + if cv_map_as_csdm: + self.cv_map_to_csdm() + + def cv_map_to_csdm(self): + # convert cv_map to csdm + self.cv_map = cp.as_csdm(np.squeeze(self.cv_map.T.copy())) + 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: + return + + d1 = cp.as_dimension(-np.log10(self.cv_lambdas), label="-log(Îģ)") + if len(self.cv_alphas) != 1: + self.cv_map.dimensions[1] = d1 + else: + self.cv_map.dimensions[0] = d1 + + def cv_plot(self): + plt.figure(figsize=(5, 3.5)) + ax = plt.subplot(projection="csdm") + ax.contour(np.log10(self.cv_map), levels=25) + ax.scatter( + -np.log10(self.hyperparameters["alpha"]), + -np.log10(self.hyperparameters["lambda"]), + marker="x", + color="k", + ) + plt.tight_layout(pad=0.5) + plt.show() + + def _get_minimizer(self): + """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=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=True, + selection="random", + ) + + if self.method == "gradient_decent": + return Lasso( + alpha=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=True, + selection="random", + ) + + if self.method == "lars": + return LassoLars( + alpha=self.cv_lambdas[0] / 2.0, + fit_intercept=False, + verbose=True, + # normalize=False, + precompute="auto", + max_iter=self.max_iterations, + eps=2.220446049250313e-16, + copy_X=True, + fit_path=False, + positive=self.positive, + jitter=None, + random_state=None, + ) + + def flip_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + y[dim0, dim1] = 0 + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def mirror_inversion(self): + csdm_obj = self.f.copy() + y = csdm_obj.y[0].components[0] + for idx, val in np.ndenumerate(y): + dim0 = idx[0] + dim1 = idx[1] + if dim1 < dim0: + if y[dim1, dim0] < 1e-6: + y[dim1, dim0] = val + new_obj = cp.CSDM( + dimensions=csdm_obj.x, dependent_variables=[cp.as_dependent_variable(y)] + ) + return new_obj + + def predict(self, K): + r"""Predict the signal using the linear model. + + Args: + K: A :math:`m \times n` kernel matrix, :math:`{\bf K}`. A numpy array of + shape (m, n). + + 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 predicted 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) + + @property + def cross_validation_curve(self): + """The cross-validation error metric determined as the mean square error. + + Returns: A two-dimensional CSDM object.""" + return self.cv_map + + +def cv(l1, X, y, cv): + """Return the cross-validation score as negative of mean square error.""" + if isinstance(l1, (Lasso, MultiTaskLasso)): + fit_params = {"check_input": False} + if isinstance(l1, LassoLars): + fit_params = None # {"Xy": np.dot(X.T, y)} + + cv_score = cross_validate( + l1, + X=X, + y=y, + scoring="neg_mean_squared_error", # 'neg_mean_absolute_error", + cv=cv, + params=fit_params, + n_jobs=1, + verbose=0, + ) + return cv_score["test_score"].mean() + + +def _get_smooth_size(f_shape, regularizer, max_size): + r"""Return the number of rows appended to for the augmented kernel. + + For smooth-lasso, the number of rows is given as + rows = \prod_{i=1}^d n_i (\sum_{j=0}^d (n_j-1)/n_j) + + For sparse ridge fusion, the number of rows is given as + rows = \prod_{i=1}^d n_i (\sum_{j=0}^d (n_j-2)/n_j) + """ + shape = np.asarray(f_shape) + shape_prod = shape.prod() + if shape_prod != max_size: + raise ValueError( + "The product of the shape must be equal to the length of axis 1 of the " + "kernel, K" + ) + if regularizer == "smooth lasso": + smooth_array_size = [int(shape_prod * (i - 1) / i) for i in shape] + smooth_size = np.asarray(smooth_array_size).sum() + elif regularizer == "sparse ridge fusion": + smooth_array_size = [int(shape_prod * (i - 2) / i) for i in shape] + smooth_size = np.asarray(smooth_array_size).sum() + else: + smooth_size = 0 + + return smooth_size + + +def _get_cv_indexes(K, folds, regularizer, f_shape=None, random=False, times=1): + """Return the indexes of the kernel and signal, corresponding to the test + and train sets.""" + cv_indexes = [] + ks0, ks1 = K.shape + if isinstance(f_shape, int): + f_shape = (f_shape,) + + tr_ = [] + smooth_size = _get_smooth_size(f_shape, regularizer, ks1) + + tr_ = (np.arange(smooth_size) + ks0).tolist() + + for j in range(folds): + train_ = [] + test_ = [] + for i in range(ks0): + if (i + j - folds) % folds == 0: + test_.append(i) + else: + train_.append(i) + train_ += tr_ + cv_indexes.append([train_, test_]) + + if random: + for _ in range(times): + kf = KFold(n_splits=folds, shuffle=True) + kf.get_n_splits(K) + for train_index, test_index in kf.split(K): + cv_indexes.append([train_index.tolist() + tr_, test_index]) + return cv_indexes + + +def generate_J_i(Ai, alpha, f_shape): + J = [] + sqrt_alpha = np.sqrt(alpha) + identity = [np.eye(i) for i in f_shape] + for i, i_count in enumerate(f_shape): + J_array = deepcopy(identity) + J_array[i] = Ai(i_count) + Ji_ = 1 + for j_array_ in J_array: + Ji_ = np.kron(Ji_, j_array_) + J.append(Ji_ * sqrt_alpha) + return J + + +def _get_augmented_data(K, s, alpha, regularizer, f_shape=None): + """Creates a smooth kernel, K, with alpha regularization parameter.""" + if alpha == 0: + return np.asfortranarray(K), np.asfortranarray(s) + + ks0, ks1 = K.shape + ss0, ss1 = s.shape + + if isinstance(f_shape, int): + f_shape = (f_shape,) + + smooth_size = _get_smooth_size(f_shape, regularizer, ks1) + + if regularizer == "smooth lasso": + + def Ai_smooth_lasso(i): + return (-1 * np.eye(i) + np.diag(np.ones(i - 1), k=-1))[1:] + + J = generate_J_i(Ai_smooth_lasso, alpha, f_shape) + + if regularizer == "sparse ridge fusion": + + def Ai_sparse_ridge_fusion(i): + A_temp = -1 * np.eye(i) + A_temp += 2 * np.diag(np.ones(i - 1), k=-1) + A_temp += -1 * np.diag(np.ones(i - 2), k=-2) + return A_temp[2:] + # return ( + # -1 * np.eye(i) + # + 2 * np.diag(np.ones(i - 1), k=-1) + # - 1 * np.diag(np.ones(i - 2), k=-2) + # )[2:] + + J = generate_J_i(Ai_sparse_ridge_fusion, alpha, f_shape) + + K_ = np.empty((ks0 + smooth_size, ks1)) + + K_[:ks0] = K + start = ks0 + end = ks0 + for J_i in J: + end = end + J_i.shape[0] + K_[start:end] = J_i + start = end + + s_ = np.zeros((ss0 + smooth_size, ss1)) + s_[:ss0] = s.real + + return np.asfortranarray(K_), np.asfortranarray(s_) + + +def prepare_signal(sig): + """Prepare signal for inversion""" + s_ = sig.y[0].components[0].T if isinstance(sig, cp.CSDM) else sig + s_ = s_[:, np.newaxis] if s_.ndim == 1 else s_ + + scale = np.sqrt(np.mean(np.abs(s_) ** 2)) + s_ = s_ / scale + return s_, scale diff --git a/mrinversion/linear_model/linear_inversion.py b/mrinversion/linear_model/linear_inversion.py index 363644f..bf8a01f 100644 --- a/mrinversion/linear_model/linear_inversion.py +++ b/mrinversion/linear_model/linear_inversion.py @@ -1,44 +1,46 @@ -import numpy as np - -__author__ = "Deepansh J. Srivastava" -__email__ = "srivastava.89@osu.edu" - - -def find_optimum_singular_value(s): - length = s.size - s2 = s**2.0 - sj = s2 / s2.sum() - T = sj * np.log10(sj) - T[np.where(np.isnan(T))] = 0 - log_n = np.log10(length) - log_nm1 = np.log10(length - 1.0) - entropy = (-1.0 / log_n) * T.sum() - - deltaEntropy = entropy - (entropy * log_n + T) / log_nm1 - - c = deltaEntropy.mean() - d = deltaEntropy.std() - - r = np.argmin(deltaEntropy - c + d) - return r - - -def TSVD(K): - U, S, VT = np.linalg.svd(K, full_matrices=False) - r = find_optimum_singular_value(S) - return U, S, VT, r - - -# standard deviation of noise remains unchanged after unitary transformation. -def reduced_subspace_kernel_and_data(U, S, VT, signal, sigma=None): - diagS = np.diag(S) - K_tilde = np.dot(diagS, VT) - s_tilde = np.dot(U.T, signal) - - projectedSignal = np.dot(U, s_tilde) - guess_solution = np.dot(np.dot(VT.T, np.diag(1 / S)), s_tilde) - - K_tilde = np.asfortranarray(K_tilde) - s_tilde = np.asfortranarray(s_tilde) - projectedSignal = np.asfortranarray(projectedSignal) - return K_tilde, s_tilde, projectedSignal, guess_solution +import numpy as np + +# import numpy.ma as ma + +__author__ = "Deepansh J. Srivastava" +__email__ = "srivastava.89@osu.edu" + + +def find_optimum_singular_value(s): + length = s.size + s2 = s**2.0 + sj = s2 / s2.sum() + T = sj * np.log10(sj) + T[np.where(np.isnan(T))] = 0 + log_n = np.log10(length) + log_nm1 = np.log10(length - 1.0) + entropy = (-1.0 / log_n) * T.sum() + + deltaEntropy = entropy - (entropy * log_n + T) / log_nm1 + + c = deltaEntropy.mean() + d = deltaEntropy.std() + + r = np.argmin(deltaEntropy - c + d) + return r + + +def TSVD(K): + U, S, VT = np.linalg.svd(K, full_matrices=False) + r = find_optimum_singular_value(S) + return U, S, VT, r + + +# standard deviation of noise remains unchanged after unitary transformation. +def reduced_subspace_kernel_and_data(U, S, VT, signal, sigma=None): + diagS = np.diag(S) + K_tilde = np.dot(diagS, VT) + s_tilde = np.dot(U.T, signal) + + projectedSignal = np.dot(U, s_tilde) + guess_solution = np.dot(np.dot(VT.T, np.diag(1 / S)), s_tilde) + + K_tilde = np.asfortranarray(K_tilde) + s_tilde = np.asfortranarray(s_tilde) + projectedSignal = np.asfortranarray(projectedSignal) + return K_tilde, s_tilde, projectedSignal, guess_solution diff --git a/mrinversion/linear_model/smooth_lasso.py b/mrinversion/linear_model/smooth_lasso.py index 9f9b202..297ff72 100644 --- a/mrinversion/linear_model/smooth_lasso.py +++ b/mrinversion/linear_model/smooth_lasso.py @@ -1,218 +1,240 @@ -from ._base_l1l2 import GeneralL2Lasso -from ._base_l1l2 import GeneralL2LassoCV - -__author__ = "Deepansh J. Srivastava" -__email__ = "srivastava.89@osu.edu" - - -class SmoothLasso(GeneralL2Lasso): - r""" - The linear model trained with the combined l1 and l2 priors as the regularizer. - The method minimizes the objective function, - - .. math:: - \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 - + \lambda \| {\bf f} \|_1 , - - where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, - :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known (measured) - signal, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` - is the desired solution. The parameters, :math:`\alpha` and :math:`\lambda`, - are the hyperparameters controlling the smoothness and sparsity of the - solution :math:`{\bf f}`. The matrix :math:`{\bf J}_i` is given as - - .. math:: - {\bf J}_i = {\bf I}_{n_1} \otimes \cdots \otimes {\bf A}_{n_i} - \otimes \cdots \otimes {\bf I}_{n_{d}}, - - where :math:`{\bf I}_{n_i} \in \mathbb{R}^{n_i \times n_i}` is the identity matrix, - - .. math:: - {\bf A}_{n_i} = \left(\begin{array}{ccccc} - 1 & -1 & 0 & \cdots & \vdots \\ - 0 & 1 & -1 & \cdots & \vdots \\ - \vdots & \vdots & \vdots & \vdots & 0 \\ - 0 & \cdots & 0 & 1 & -1 - \end{array}\right) \in \mathbb{R}^{(n_i-1)\times n_i}, - - and the symbol :math:`\otimes` is the Kronecker product. The terms, - :math:`\left(n_1, n_2, \cdots, n_d\right)`, are the number of points along the - respective dimensions, with the constraint that :math:`\prod_{i=1}^{d}n_i = n`, - where :math:`d` is the total number of dimensions. - - Args - ---- - - alpha: float - The hyperparameter, :math:`\alpha`. - lambda1: float - The hyperparameter, :math:`\lambda`. - inverse_dimension: list - A list of csdmpy Dimension objects representing the inverse space. - max_iterations: int - The maximum number of iterations allowed when solving the problem. The default - value is 10000. - tolerance: float - The tolerance at which the solution is considered converged. The default value - is 1e-5. - positive: bool - If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only - positive values, else the solution may contain positive and negative amplitudes. - The default is True. - - Attributes - ---------- - - f: ndarray or CSDM object. - A ndarray of shape (`m_count`, `nd`, ..., `n1`, `n0`) representing the - solution, :math:`{\bf f} \in \mathbb{R}^{m_\text{count} \times n_d \times - \cdots n_1 \times n_0}`. - n_iter: int - 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, - ) - - -class SmoothLassoCV(GeneralL2LassoCV): - r""" - The linear model trained with the combined l1 and l2 priors as the - regularizer. The method minimizes the objective function, - - .. math:: - \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 - + \lambda \| {\bf f} \|_1 , - - where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, - :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known signal - containing noise, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` - is the desired solution. The parameters, :math:`\alpha` and :math:`\lambda`, - are the hyperparameters controlling the smoothness and sparsity of the - solution :math:`{\bf f}`. - The matrix :math:`{\bf J}_i` is given as - - .. math:: - {\bf J}_i = {\bf I}_{n_1} \otimes \cdots \otimes {\bf A}_{n_i} - \otimes \cdots \otimes {\bf I}_{n_{d}}, - - where :math:`{\bf I}_{n_i} \in \mathbb{R}^{n_i \times n_i}` is the identity - matrix, - - .. math:: - {\bf A}_{n_i} = \left(\begin{array}{ccccc} - 1 & -1 & 0 & \cdots & \vdots \\ - 0 & 1 & -1 & \cdots & \vdots \\ - \vdots & \vdots & \vdots & \vdots & 0 \\ - 0 & \cdots & 0 & 1 & -1 - \end{array}\right) \in \mathbb{R}^{(n_i-1)\times n_i}, - - and the symbol :math:`\otimes` is the Kronecker product. The terms, - :math:`\left(n_1, n_2, \cdots, n_d\right)`, are the number of points along the - 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. - inverse_dimension: list - A list of csdmpy Dimension objects representing the inverse space. - folds: int - The number of folds used in cross-validation.The default is 10. - max_iterations: int - The maximum number of iterations allowed when solving the problem. The default - value is 10000. - tolerance: float - The tolerance at which the solution is considered converged. The default value - is 1e-5. - positive: bool - If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only - positive values, else the solution may contain positive and negative amplitudes. - The default is True. - 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. - times: int - The number of times to randomized n-folds are created. Only applicable when - `randomize` attribute is True. - verbose: bool - If true, prints the process. - n_jobs: int - The number of CPUs used for computation. The default is -1, that is, all - available CPUs are used. - - - Attributes - ---------- - f: ndarray or CSDM object. - A ndarray of shape (m_count, nd, ..., n1, n0). The solution, - :math:`{\bf f} \in \mathbb{R}^{m_\text{count} \times n_d \times \cdots n_1 - \times n_0}` or an equivalent CSDM object. - n_iter: int. - The number of iterations required to reach the specified tolerance. - hyperparameters: dict. - A dictionary with the :math:`\alpha` and :math:\lambda` hyperparameters. - cross_validation_curve: CSDM object. - 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, - ) +from ._base_l1l2 import GeneralL2Lasso +from ._base_l1l2 import GeneralL2LassoCV + +__author__ = "Deepansh J. Srivastava" +__email__ = "srivastava.89@osu.edu" + + +class SmoothLasso(GeneralL2Lasso): + r""" + The linear model trained with the combined l1 and l2 priors as the regularizer. + The method minimizes the objective function, + + .. math:: + \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 + + \lambda \| {\bf f} \|_1 , + + where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, + :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known (measured) + signal, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` + is the desired solution. The parameters, :math:`\alpha` and :math:`\lambda`, + are the hyperparameters controlling the smoothness and sparsity of the + solution :math:`{\bf f}`. The matrix :math:`{\bf J}_i` is given as + + .. math:: + {\bf J}_i = {\bf I}_{n_1} \otimes \cdots \otimes {\bf A}_{n_i} + \otimes \cdots \otimes {\bf I}_{n_{d}}, + + where :math:`{\bf I}_{n_i} \in \mathbb{R}^{n_i \times n_i}` is the identity matrix, + + .. math:: + {\bf A}_{n_i} = \left(\begin{array}{ccccc} + 1 & -1 & 0 & \cdots & \vdots \\ + 0 & 1 & -1 & \cdots & \vdots \\ + \vdots & \vdots & \vdots & \vdots & 0 \\ + 0 & \cdots & 0 & 1 & -1 + \end{array}\right) \in \mathbb{R}^{(n_i-1)\times n_i}, + + and the symbol :math:`\otimes` is the Kronecker product. The terms, + :math:`\left(n_1, n_2, \cdots, n_d\right)`, are the number of points along the + respective dimensions, with the constraint that :math:`\prod_{i=1}^{d}n_i = n`, + where :math:`d` is the total number of dimensions. + + Args + ---- + + alpha: float + The hyperparameter, :math:`\alpha`. + lambda1: float + The hyperparameter, :math:`\lambda`. + inverse_dimension: list + A list of csdmpy Dimension objects representing the inverse space. + max_iterations: int + The maximum number of iterations allowed when solving the problem. The default + value is 10000. + tolerance: float + The tolerance at which the solution is considered converged. The default value + is 1e-5. + positive: bool + If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only + positive values, else the solution may contain positive and negative amplitudes. + The default is True. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. + + Attributes + ---------- + + f: ndarray or CSDM object. + A ndarray of shape (`m_count`, `nd`, ..., `n1`, `n0`) representing the + solution, :math:`{\bf f} \in \mathbb{R}^{m_\text{count} \times n_d \times + \cdots n_1 \times n_0}`. + n_iter: int + 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", + xygrid="full", + ): + super().__init__( + alpha=alpha, + lambda1=lambda1, + max_iterations=max_iterations, + tolerance=tolerance, + positive=positive, + regularizer="smooth lasso", + inverse_dimension=inverse_dimension, + method=method, + xygrid=xygrid, + ) + + +class SmoothLassoCV(GeneralL2LassoCV): + r""" + The linear model trained with the combined l1 and l2 priors as the + regularizer. The method minimizes the objective function, + + .. math:: + \| {\bf Kf - s} \|^2_2 + \alpha \sum_{i=1}^{d} \| {\bf J}_i {\bf f} \|_2^2 + + \lambda \| {\bf f} \|_1 , + + where :math:`{\bf K} \in \mathbb{R}^{m \times n}` is the kernel, + :math:`{\bf s} \in \mathbb{R}^{m \times m_\text{count}}` is the known signal + containing noise, and :math:`{\bf f} \in \mathbb{R}^{n \times m_\text{count}}` + is the desired solution. The parameters, :math:`\alpha` and :math:`\lambda`, + are the hyperparameters controlling the smoothness and sparsity of the + solution :math:`{\bf f}`. + The matrix :math:`{\bf J}_i` is given as + + .. math:: + {\bf J}_i = {\bf I}_{n_1} \otimes \cdots \otimes {\bf A}_{n_i} + \otimes \cdots \otimes {\bf I}_{n_{d}}, + + where :math:`{\bf I}_{n_i} \in \mathbb{R}^{n_i \times n_i}` is the identity + matrix, + + .. math:: + {\bf A}_{n_i} = \left(\begin{array}{ccccc} + 1 & -1 & 0 & \cdots & \vdots \\ + 0 & 1 & -1 & \cdots & \vdots \\ + \vdots & \vdots & \vdots & \vdots & 0 \\ + 0 & \cdots & 0 & 1 & -1 + \end{array}\right) \in \mathbb{R}^{(n_i-1)\times n_i}, + + and the symbol :math:`\otimes` is the Kronecker product. The terms, + :math:`\left(n_1, n_2, \cdots, n_d\right)`, are the number of points along the + 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. + inverse_dimension: list + A list of csdmpy Dimension objects representing the inverse space. + folds: int + The number of folds used in cross-validation.The default is 10. + max_iterations: int + The maximum number of iterations allowed when solving the problem. The default + value is 10000. + tolerance: float + The tolerance at which the solution is considered converged. The default value + is 1e-5. + positive: bool + If True, the amplitudes in the solution, :math:`{\bf f}`, is constrained to only + positive values, else the solution may contain positive and negative amplitudes. + The default is True. + 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. + times: int + The number of times to randomized n-folds are created. Only applicable when + `randomize` attribute is True. + verbose: bool + If true, prints the process. + n_jobs: int + The number of CPUs used for computation. The default is -1, that is, all + available CPUs are used. + xygrid: 'full', 'mirrored', 'positive', or 'negative' + If 'full', the inversion is performed over the entire xygrid and returned as-is. + For quadrupolar kernels, 'mirrored', 'positive', or 'negative' should be used. + If 'mirrored', the inversion is calculated over the positive half-quadrant and + reflected over the :math:`\eta_q=1` line to display both :math:`+C_q` and + :math:`-C_q` results. If 'positive', the inversion is calculated over the + positive half-quadrant and returned as is. If 'negative', the inversion is + calculated over the positive half-quadrant and flipped about the + :math:`\eta_q=1` line to display :math:`-C_q` results. The default is 'full'. + + + Attributes + ---------- + f: ndarray or CSDM object. + A ndarray of shape (m_count, nd, ..., n1, n0). The solution, + :math:`{\bf f} \in \mathbb{R}^{m_\text{count} \times n_d \times \cdots n_1 + \times n_0}` or an equivalent CSDM object. + n_iter: int. + The number of iterations required to reach the specified tolerance. + hyperparameters: dict. + A dictionary with the :math:`\alpha` and :math:\lambda` hyperparameters. + cross_validation_curve: CSDM object. + 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", + xygrid="full", + ): + 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, + xygrid=xygrid, + ) diff --git a/mrinversion/linear_model/utils.py b/mrinversion/linear_model/utils.py new file mode 100644 index 0000000..a65d071 --- /dev/null +++ b/mrinversion/linear_model/utils.py @@ -0,0 +1,226 @@ +import numpy as np +from scipy.optimize import minimize + +from mrinversion.linear_model import SmoothLassoCV + + +class CVMinimizer: + def __init__( + self, + inverse_dimension, + compressed_K, + compressed_s, + guess_neglogalpha, + guess_negloglambda, + sigma=0.0, + folds=10, + xygrid="full", + ): + self.inverse_dimension = inverse_dimension + self.compressed_K = compressed_K + self.compressed_s = compressed_s + self.guess = np.array([guess_neglogalpha, guess_negloglambda]) + self.sigma = sigma + self.folds = folds + self.xygrid = xygrid + + @staticmethod + def _cv_min_withprint(loghyperparameters, self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid, + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + print(f"alpha: {-np.log10(alpha)}, lambda: {-np.log10(lam)}, cv: {cv}") + return cv + + @staticmethod + def _cv_min_noprint(loghyperparameters, self): + alpha = 10 ** (-loghyperparameters[0]) + lam = 10 ** (-loghyperparameters[1]) + + s_lasso_fit = SmoothLassoCV( + alphas=np.array([alpha]), + lambdas=np.array([lam]), + inverse_dimension=self.inverse_dimension, + sigma=self.sigma, + folds=self.folds, + xygrid=self.xygrid, + ) + + s_lasso_fit.fit(self.compressed_K, self.compressed_s, cv_map_as_csdm=False) + cv = np.log10(s_lasso_fit.cv_map[0][0]) + return cv + + @staticmethod + def _print_result(this_result, print_steps, prev_result=None): + if prev_result: + if this_result.fun < prev_result.fun: + if this_result.success: + print( + "This minimization improved on the previous result and \n" + "completed the minimization in the iterations given. " + ) + else: + print( + "This minimization improved on the previous result, but did \n" + "not finish in the given number of iterations." + ) + else: + if this_result.success: + print("This minimization minimized to a worse solution.") + else: + print( + "This minimization did not improve on the previous result, \n" + "and did not finish in the given number of iterations." + ) + if print_steps: + print(f"This minimization results: \n\n{this_result}") + else: + if this_result.success: + print("This minimization completed in the iterations given. ") + else: + print("This minimization did not finish in given number of iterations.") + if print_steps: + print(f"This minimization results: \n\n{this_result}") + + print("--------------------------------------------------\n") + + def _minimize( + self, guess, tol, maxiter, print_funcevals, print_steps, prev_result=None + ): + if print_funcevals: + simplex_res = minimize( + fun=self._cv_min_withprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={"maxiter": maxiter}, + method="Nelder-Mead", + ) + else: + simplex_res = minimize( + fun=self._cv_min_noprint, + x0=guess, + args=(self), + tol=tol, + # bounds=((2,10), (2,10)), + options={"maxiter": maxiter}, + method="Nelder-Mead", + ) + self._print_result( + this_result=simplex_res, print_steps=print_steps, prev_result=prev_result + ) + return simplex_res + + def mincv_close( + self, tol=0.0002, maxiter=500, print_steps=True, print_funcevals=False + ): + simplex_res1 = self._minimize( + guess=self.guess, + tol=tol * 1.5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + ) + + if simplex_res1.success: + simplex_res2 = self._minimize( + guess=simplex_res1.x, + tol=tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + prev_result=simplex_res1, + ) + if simplex_res2.fun < simplex_res1.fun: + return simplex_res2 + else: + return simplex_res1 + else: + return simplex_res1 + + def mincv( + self, + tol=0.0002, + maxiter=100, + print_steps=True, + print_funcevals=False, + guess_bounds=[[2, 10], [2, 10]], + guesses="all", + ): + '''options for guesses: "all", "cartesian", "diagonals"''' + lowval_alpha = np.average([self.guess[0], guess_bounds[0][0]]) + highval_alpha = np.average([self.guess[0], guess_bounds[0][1]]) + lowval_lambda = np.average([self.guess[1], guess_bounds[1][0]]) + highval_lambda = np.average([self.guess[1], guess_bounds[1][1]]) + if guesses == "all": + trythese = np.zeros((9, 2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, highval_lambda] + trythese[2] = [highval_alpha, self.guess[1]] + trythese[3] = [highval_alpha, lowval_lambda] + trythese[4] = [self.guess[0], lowval_lambda] + trythese[5] = [lowval_alpha, lowval_lambda] + trythese[6] = [lowval_alpha, self.guess[1]] + trythese[7] = [lowval_alpha, highval_lambda] + trythese[8] = self.guess + elif guesses == "cartesian": + trythese = np.zeros((5, 2)) + trythese[0] = [self.guess[0], highval_lambda] + trythese[1] = [highval_alpha, self.guess[1]] + trythese[2] = [self.guess[0], lowval_lambda] + trythese[3] = [lowval_alpha, self.guess[1]] + trythese[4] = self.guess + elif guesses == "diagonal": + trythese = np.zeros((5, 2)) + trythese[0] = [highval_alpha, highval_lambda] + trythese[1] = [highval_alpha, lowval_lambda] + trythese[2] = [lowval_alpha, lowval_lambda] + trythese[3] = [lowval_alpha, highval_lambda] + trythese[4] = self.guess + else: + print('choose a valid choice for "guesses"') + return None + print(f"Starting points to try: {trythese}") + first_step = [ + self._minimize( + guess=thisguess, + tol=tol * 5, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=print_steps, + ) + for thisguess in trythese + ] + print("-----------------------------------------------------------------") + mins = np.asarray([step.fun for step in first_step]) + min_idx = np.where(mins == mins.min())[0][0] + + final_step = self._minimize( + guess=first_step[min_idx].x, + tol=tol, + maxiter=maxiter, + print_funcevals=print_funcevals, + print_steps=False, + ) + + self._print_result( + final_step, print_steps=True, prev_result=first_step[min_idx] + ) + + if final_step.fun > first_step[min_idx].fun: + print("Second minimization did not improve.") + return first_step[min_idx] + else: + return final_step diff --git a/mrinversion/tests/fista_intergration_test.py b/mrinversion/tests/fista_intergration_test.py index 03dbe5d..1290756 100644 --- a/mrinversion/tests/fista_intergration_test.py +++ b/mrinversion/tests/fista_intergration_test.py @@ -1,55 +1,55 @@ -"""Integration tests for FISTA implementation.""" -import csdmpy as cp -import numpy as np - -from mrinversion.kernel import relaxation -from mrinversion.linear_model import LassoFistaCV -from mrinversion.linear_model import TSVDCompression - - -def test_fista(): - domain = "https://ssnmr.org/resources/mrinversion" - filename = f"{domain}/test1_signal.csdf" - signal = cp.load(filename) - sigma = 0.0008 - - datafile = f"{domain}/test1_t2.csdf" - true_dist = cp.load(datafile) - - kernel_dimension = signal.dimensions[0] - relaxT2 = relaxation.T2( - kernel_dimension=kernel_dimension, - inverse_dimension=dict( - count=64, - minimum="1e-2 s", - maximum="1e3 s", - scale="log", - label="log (T2 / s)", - ), - ) - inverse_dimension = relaxT2.inverse_dimension - K = relaxT2.kernel(supersampling=1) - - new_system = TSVDCompression(K, signal) - compressed_K = new_system.compressed_K - compressed_s = new_system.compressed_s - - assert new_system.truncation_index == 29 - - lambdas = 10 ** (-5 + 4 * (np.arange(16) / 15)) - - f_lasso_cv = LassoFistaCV( - lambdas=lambdas, # A numpy array of lambda values. - folds=5, # The number of folds in n-folds cross-validation. - sigma=sigma, # noise standard deviation - inverse_dimension=inverse_dimension, # previously defined inverse dimensions. - ) - - f_lasso_cv.fit(K=compressed_K, s=compressed_s) - - sol = f_lasso_cv.f - assert np.argmax(sol.y[0].components[0]) == np.argmax(true_dist.y[0].components[0]) - - residuals = f_lasso_cv.residuals(K=K, s=signal) - std = residuals.std() - np.testing.assert_almost_equal(std.value, sigma, decimal=2) +"""Integration tests for FISTA implementation.""" +import csdmpy as cp +import numpy as np + +from mrinversion.kernel import relaxation +from mrinversion.linear_model import LassoFistaCV +from mrinversion.linear_model import TSVDCompression + + +def test_fista(): + domain = "https://ssnmr.org/resources/mrinversion" + filename = f"{domain}/test1_signal.csdf" + signal = cp.load(filename) + sigma = 0.0008 + + datafile = f"{domain}/test1_t2.csdf" + true_dist = cp.load(datafile) + + kernel_dimension = signal.dimensions[0] + relaxT2 = relaxation.T2( + kernel_dimension=kernel_dimension, + inverse_dimension=dict( + count=64, + minimum="1e-2 s", + maximum="1e3 s", + scale="log", + label="log (T2 / s)", + ), + ) + inverse_dimension = relaxT2.inverse_dimension + K = relaxT2.kernel(supersampling=1) + + new_system = TSVDCompression(K, signal) + compressed_K = new_system.compressed_K + compressed_s = new_system.compressed_s + + assert new_system.truncation_index == 29 + + lambdas = 10 ** (-5 + 4 * (np.arange(16) / 15)) + + f_lasso_cv = LassoFistaCV( + lambdas=lambdas, # A numpy array of lambda values. + folds=5, # The number of folds in n-folds cross-validation. + sigma=sigma, # noise standard deviation + inverse_dimension=inverse_dimension, # previously defined inverse dimensions. + ) + + f_lasso_cv.fit(K=compressed_K, s=compressed_s) + + sol = f_lasso_cv.f + assert np.argmax(sol.y[0].components[0]) == np.argmax(true_dist.y[0].components[0]) + + residuals = f_lasso_cv.residuals(K=K, s=signal) + std = residuals.std() + np.testing.assert_almost_equal(std.value, sigma, decimal=2) diff --git a/mrinversion/utils.py b/mrinversion/utils.py index c4b55b9..79b8c52 100644 --- a/mrinversion/utils.py +++ b/mrinversion/utils.py @@ -1,451 +1,623 @@ -from copy import deepcopy -from itertools import combinations -from itertools import product - -import csdmpy as cp -import matplotlib.pyplot as plt -import numpy as np -from matplotlib import cm -from mpl_toolkits.mplot3d import Axes3D # lgtm [py/unused-import] - - -def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): - """Convert the three-dimensional p(iso, x, y) to p(iso, zeta, eta) tensor - distribution. - - Args - ---- - - csdm_object: CSDM - A CSDM object containing the 3D p(iso, x, y) distribution. - zeta: CSDM.Dimension - A CSDM dimension object describing the zeta dimension. - eta: CSDM.Dimension - A CSDM dimension object describing the eta dimension. - n: int - An integer used in linear interpolation of the data. The default is 5. - """ - [ - item.to("ppm", "nmr_frequency_ratio") - for item in csdm_object.x - if item.origin_offset != 0 - ] - data = csdm_object.y[0].components[0] - - extra_dims = 1 - if len(csdm_object.x) > 2: - extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) - data.shape = (extra_dims, data.shape[-2], data.shape[-1]) - - reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) - dx = reg_x[1] - reg_x[0] - dy = reg_y[1] - reg_y[0] - sol = np.zeros((extra_dims, zeta.count, eta.count)) - - bins = [zeta.count, eta.count] - d_zeta = zeta.increment.value / 2 - d_eta = eta.increment.value / 2 - range_ = [ - [zeta.coordinates[0].value - d_zeta, zeta.coordinates[-1].value + d_zeta], - [eta.coordinates[0] - d_eta, eta.coordinates[-1] + d_eta], - ] - - avg_range_x = (np.arange(n) - (n - 1) / 2) * dx / n - avg_range_y = (np.arange(n) - (n - 1) / 2) * dy / n - for x_item in avg_range_x: - for y_item in avg_range_y: - x__ = np.abs(reg_x + x_item) - y__ = np.abs(reg_y + y_item) - x_, y_ = np.meshgrid(x__, y__) - x_ = x_.ravel() - y_ = y_.ravel() - - zeta_grid = np.sqrt(x_**2 + y_**2) - eta_grid = np.ones(zeta_grid.shape) - - index = np.where(x_ < y_) - eta_grid[index] = (4 / np.pi) * np.arctan(x_[index] / y_[index]) - - index = np.where(x_ > y_) - zeta_grid[index] *= -1 - eta_grid[index] = (4 / np.pi) * np.arctan(y_[index] / x_[index]) - - index = np.where(x_ == y_) - np.append(zeta, -zeta_grid[index]) - np.append(eta, np.ones(index[0].size)) - for i in range(extra_dims): - weight = deepcopy(data[i]).ravel() - weight[index] /= 2 - np.append(weight, weight[index]) - sol_, _, _ = np.histogram2d( - zeta_grid, eta_grid, weights=weight, bins=bins, range=range_ - ) - sol[i] += sol_ - - sol /= n * n - - del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y - csdm_new = cp.as_csdm(np.squeeze(sol)) - csdm_new.x[0] = eta - csdm_new.x[1] = zeta - if len(csdm_object.x) > 2: - csdm_new.x[2] = csdm_object.x[2] - return csdm_new - - -def get_polar_grids(ax, ticks=None, offset=0): - """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. - - Args: - ax: Matplotlib Axes. - ticks: Tick coordinates where radial grids are drawn. The value can be a list - or a numpy array. The default value is None. - offset: The grid is drawn at an offset away from the origin. - """ - ylim = ax.get_ylim() - xlim = ax.get_xlim() - if ticks is None: - x = np.asarray(ax.get_xticks()) - inc = x[1] - x[0] - size = x.size - x = np.arange(size + 5) * inc - else: - x = np.asarray(ticks) - - lw = 0.3 - t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) - t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) - - ax.add_artist(t1) - ax.add_artist(t2) - for x_ in x: - if x_ - offset > 0: - ax.add_artist( - plt.Circle( - (0, 0), - x_ - offset, - fill=False, - color="k", - linestyle="--", - linewidth=lw, - alpha=0.5, - ) - ) - - angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 4.0) - angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) - for ang_ in angle1: - ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) - for ang_ in angle2: - ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) - ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) - ax.set_xlim(xlim) - ax.set_ylim(ylim) - - -def plot_3d( - ax, - csdm_objects, - elev=28, - azim=-150, - x_lim=None, - y_lim=None, - z_lim=None, - cmap=cm.PiYG, - box=False, - clip_percent=0.0, - linewidth=1, - alpha=0.15, - **kwargs, -): - r"""Generate a 3D density plot with 2D contour and 1D projections. - - Args: - ax: Matplotlib Axes to render the plot. - csdm_objects: A 3D{1} CSDM object or a list of CSDM objects holding the data. - elev: (optional) The 3D view angle, elevation angle in the z plane. - azim: (optional) The 3D view angle, azimuth angle in the x-y plane. - x_lim: (optional) The x limit given as a list, [x_min, x_max]. - y_lim: (optional) The y limit given as a list, [y_min, y_max]. - z_lim: (optional) The z limit given as a list, [z_min, z_max]. - max_2d: (Optional) The normalization factor of the 2D contour projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the - maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. - max_1d: (Optional) The normalization factor of the 1D projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the - maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. - cmap: (Optional) The colormap or a list of colormaps used in rendering the - volumetric plot. The colormap list is applied to the ordered list of csdm_objects. - The same colormap is used for the 2D contour projections. For 1D plots, the first - color in the colormap scheme is used for the line color. - box: (Optional) If True, draw a box around the 3D data region. - clip_percent: (Optional) The amplitudes of the dataset below the given percent - is made transparent for the volumetric plot. - linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. - alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D - volume. - """ - csdm_object_list = csdm_objects - if not isinstance(csdm_objects, list): - csdm_object_list = [csdm_objects] - - if not isinstance(cmap, list): - cmap = [cmap] - - fraction = np.array([abs(item).sum() for item in csdm_object_list]) - fraction /= fraction.max() - - for index, item in enumerate(csdm_object_list): - plot_3d_( - ax, - item, - elev=elev, - azim=azim, - x_lim=x_lim, - y_lim=y_lim, - z_lim=z_lim, - fraction=fraction[index], - cmap=cmap[index], - box=box, - clip_percent=clip_percent, - linewidth=linewidth, - alpha=alpha, - **kwargs, - ) - - -def plot_3d_( - ax, - csdm_object, - elev=28, - azim=-150, - x_lim=None, - y_lim=None, - z_lim=None, - fraction=1.0, - cmap=cm.PiYG, - box=False, - clip_percent=0.0, - linewidth=1, - alpha=0.15, - **kwargs, -): - r"""Generate a 3D density plot with 2D contour and 1D projections. - - Args: - ax: Matplotlib Axes to render the plot. - csdm_object: A 3D{1} CSDM object holding the data. - elev: (optional) The 3D view angle, elevation angle in the z plane. - azim: (optional) The 3D view angle, azimuth angle in the x-y plane. - x_lim: (optional) The x limit given as a list, [x_min, x_max]. - y_lim: (optional) The y limit given as a list, [y_min, y_max]. - z_lim: (optional) The z limit given as a list, [z_min, z_max]. - max_2d: (Optional) The normalization factor of the 2D contour projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the - maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. - max_1d: (Optional) The normalization factor of the 1D projections. The - attribute is meaningful when multiple 3D datasets are viewed on the same - plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the - maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. - cmap: (Optional) The colormap used in rendering the volumetric plot. The same - colormap is used for the 2D contour projections. For 1D plots, the first - color in the colormap scheme is used for the line color. - box: (Optional) If True, draw a box around the 3D data region. - clip_percent: (Optional) The amplitudes of the dataset below the given percent - is made transparent for the volumetric plot. - linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. - alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D - volume. - """ - lw = linewidth - if isinstance(csdm_object, cp.CSDM): - f = csdm_object.y[0].components[0].T - - a_, b_, c_ = (item for item in csdm_object.x) - - a = a_.coordinates.value - b = b_.coordinates.value - c = c_.coordinates.value - - xlabel = f"{a_.axis_label} - 0" - ylabel = f"{b_.axis_label} - 1" - zlabel = f"{c_.axis_label} - 2" - - else: - f = csdm_object - a = np.arange(f.shape[0]) - b = np.arange(f.shape[1]) - c = np.arange(f.shape[2]) - - xlabel = "x" - ylabel = "y" - zlabel = "z" - - delta_a = a[1] - a[0] - delta_b = b[1] - b[0] - delta_c = c[1] - c[0] - - clr = cmap - ck = cmap(0) - facecolors = cmap(f) - facecolors[:, :, :, -1] = f * alpha - index = np.where(f < clip_percent / 100) - facecolors[:, :, :, -1][index] = 0 - facecolors.shape = (f.size, 4) - - if x_lim is None: - x_lim = [a[0], a[-1]] - if y_lim is None: - y_lim = [b[0], b[-1]] - if z_lim is None: - z_lim = [c[0], c[-1]] - - offset_scalar = 50 - height_scalar = 10 - z_offset = (z_lim[1] - z_lim[0]) / offset_scalar - z_scale = np.abs(z_lim[1] - z_lim[0]) / height_scalar - offz_n = z_lim[0] - (delta_c / 2.0) - offz_1d = z_lim[1] + z_offset - - y_offset = (y_lim[1] - y_lim[0]) / offset_scalar - y_scale = np.abs(y_lim[1] - y_lim[0]) / height_scalar - offy = y_lim[1] + (delta_b / 2.0) - offy_n_1d = y_lim[0] - y_offset - - offx = x_lim[1] + (delta_a / 2.0) - x_scale = np.abs(x_lim[1] - x_lim[0]) / height_scalar - - if azim > -90 and azim <= 0: - offx = x_lim[0] - (delta_a / 2.0) - offy_n_1d = y_lim[0] - y_offset - - if azim > 0 and azim <= 90: - offy = y_lim[0] - (delta_b / 2.0) - offy_n_1d = y_lim[1] + y_offset - offx = x_lim[0] - (delta_a / 2.0) - - if azim > 90: - offx = x_lim[1] + (delta_a / 2.0) - offy_n_1d = y_lim[1] + y_offset - offy = y_lim[0] - (delta_b / 2.0) - - ax.set_proj_type("persp") - ax.view_init(elev=elev, azim=azim) - - # 2D x-y contour projection --------------------- - levels = (np.arange(20) + 1) / 20 - x1, y1 = np.meshgrid(a, b, indexing="ij") - dist_xy = f.mean(axis=2) - dist_xy *= fraction * z_scale / np.abs(dist_xy.max()) - ax.contour( - x1, - y1, - dist_xy, - zdir="z", - offset=offz_n, - cmap=clr, - levels=z_scale * levels, - linewidths=lw, - ) - - # 2D x-z contour projection - x1, y1 = np.meshgrid(a, c, indexing="ij") - dist_xz = f.mean(axis=1) - dist_xz *= fraction * y_scale / np.abs(dist_xz.max()) - ax.contour( - x1, - dist_xz, - y1, - zdir="y", - offset=offy, - cmap=clr, - levels=y_scale * levels, - linewidths=lw, - ) - - # 1D x-axis projection from 2D x-z projection - proj_x = dist_xz.mean(axis=1) - proj_x *= fraction * z_scale / np.abs(proj_x.max()) - ax.plot(a, np.sign(z_offset) * proj_x + offz_1d, offy, zdir="y", c=ck, linewidth=lw) - - # 1D z-axis projection from 2D x-z projection - proj_z = dist_xz.mean(axis=0) - proj_z *= fraction * y_scale / np.abs(proj_z.max()) - ax.plot( - np.sign(azim) * np.sign(y_offset) * proj_z + offy_n_1d, - c, - offx, - zdir="x", - c=ck, - linewidth=lw, - ) - ax.set_xlim(z_lim) - - # 2D y-z contour projection - x1, y1 = np.meshgrid(b, c, indexing="ij") - dist_yz = f.mean(axis=0) - dist_yz *= fraction * x_scale / np.abs(dist_yz.max()) - ax.contour( - dist_yz, - x1, - y1, - zdir="x", - offset=offx, - cmap=clr, - levels=x_scale * levels, - linewidths=lw, - ) - - # 1D y-axis projection - proj_y = dist_yz.mean(axis=1) - proj_y *= fraction * z_scale / np.abs(proj_y.max()) - ax.plot(b, np.sign(z_offset) * proj_y + offz_1d, offx, zdir="x", c=ck, linewidth=lw) - - ax.set_xlim(x_lim) - ax.set_ylim(y_lim) - ax.set_zlim(z_lim) - - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.set_zlabel(zlabel) - - x, y, z = np.meshgrid(a, b, c, indexing="ij") - - if "s" not in kwargs: - kwargs["s"] = 300 - ax.scatter(x.flat, y.flat, z.flat, marker="X", c=facecolors, **kwargs) - - # full box - r1 = [x_lim[0] - delta_a / 2, x_lim[-1] + delta_a / 2] - r2 = [y_lim[0] - delta_b / 2, y_lim[-1] + delta_b / 2] - r3 = [z_lim[-1] - delta_c / 2, z_lim[0] + delta_c / 2] - - l_box = lw - for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): - if np.sum(np.abs(s - e)) == r1[1] - r1[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - if np.sum(np.abs(s - e)) == r2[1] - r2[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - if np.sum(np.abs(s - e)) == r3[1] - r3[0]: - ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) - - # draw cube - if box: - r1 = [a[0] - delta_a / 2, a[-1] + delta_a / 2] - r2 = [b[0] - delta_b / 2, b[-1] + delta_b / 2] - r3 = [c[0] - delta_c / 2, c[-1] + delta_c / 2] - - for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): - if np.sum(np.abs(s - e)) == r1[1] - r1[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - if np.sum(np.abs(s - e)) == r2[1] - r2[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - if np.sum(np.abs(s - e)) == r3[1] - r3[0]: - ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) - - ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) - ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) +from copy import deepcopy +from itertools import combinations +from itertools import product + +import csdmpy as cp +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import cm +from mpl_toolkits.mplot3d import Axes3D # lgtm [py/unused-import] + + +def to_Haeberlen_grid(csdm_object, zeta, eta, n=5): + """Convert the three-dimensional p(iso, x, y) to p(iso, zeta, eta) tensor + distribution. + + Args + ---- + + csdm_object: CSDM + A CSDM object containing the 3D p(iso, x, y) distribution. + zeta: CSDM.Dimension + A CSDM dimension object describing the zeta dimension. + eta: CSDM.Dimension + A CSDM dimension object describing the eta dimension. + n: int + An integer used in linear interpolation of the data. The default is 5. + """ + [ + item.to("ppm", "nmr_frequency_ratio") + for item in csdm_object.x + if item.origin_offset != 0 + ] + data = csdm_object.y[0].components[0] + # print(f'data max: {data.max()}') + extra_dims = 1 + if len(csdm_object.x) > 2: + extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) + data.shape = (extra_dims, data.shape[-2], data.shape[-1]) + + reg_x, reg_y = (csdm_object.x[i].coordinates.value for i in range(2)) + dx = reg_x[1] - reg_x[0] + dy = reg_y[1] - reg_y[0] + sol = np.zeros((extra_dims, zeta.count, eta.count)) + + bins = [zeta.count, eta.count] + d_zeta = zeta.increment.value / 2 + d_eta = eta.increment.value / 2 + range_ = [ + [zeta.coordinates[0].value - d_zeta, zeta.coordinates[-1].value + d_zeta], + [eta.coordinates[0].value - d_eta, eta.coordinates[-1].value + d_eta], + ] + + avg_range_x = (np.arange(n) - (n - 1) / 2) * dx / n + avg_range_y = (np.arange(n) - (n - 1) / 2) * dy / n + for x_item in avg_range_x: + for y_item in avg_range_y: + x__ = np.abs(reg_x + x_item) + y__ = np.abs(reg_y + y_item) + x_, y_ = np.meshgrid(x__, y__) + x_ = x_.ravel() + y_ = y_.ravel() + + zeta_grid = np.sqrt(x_**2 + y_**2) + eta_grid = np.ones(zeta_grid.shape) + + index = np.where(x_ < y_) + eta_grid[index] = (4 / np.pi) * np.arctan(x_[index] / y_[index]) + + index = np.where(x_ > y_) + zeta_grid[index] *= -1 + eta_grid[index] = (4 / np.pi) * np.arctan(y_[index] / x_[index]) + + index = np.where(x_ == y_) + np.append(zeta, -zeta_grid[index]) + np.append(eta, np.ones(index[0].size)) + for i in range(extra_dims): + weight = deepcopy(data[i]).ravel() + # print(f'weight max: {weight.max()}') + # print(f'range: {range_}') + weight[index] /= 2 + np.append(weight, weight[index]) + sol_, _, _ = np.histogram2d( + zeta_grid, eta_grid, weights=weight, bins=bins, range=range_ + ) + sol[i] += sol_ + # print(f'sol max: {sol.max()}') + # print() + sol /= n * n + + del zeta_grid, eta_grid, index, x_, y_, avg_range_x, avg_range_y + csdm_new = cp.as_csdm(np.squeeze(sol)) + csdm_new.x[0] = eta + csdm_new.x[1] = zeta + if len(csdm_object.x) > 2: + csdm_new.x[2] = csdm_object.x[2] + return csdm_new + + +def to_old_xq_yq_grid(csdm_object, xq, yq, n=5): + """Convert the three-dimensional p(iso, Cq, eta) to p(iso, xq, yq) tensor + distribution. + + Args + ---- + + csdm_object: CSDM + A CSDM object containing the 3D p(iso, Cq, eta) distribution. + xq: CSDM.Dimension + A CSDM dimension object describing the xq dimension. + yq: CSDM.Dimension + A CSDM dimension object describing the yq dimension. + n: int + An integer used in linear interpolation of the data. The default is 5. + """ + [ + item.to("ppm", "nmr_frequency_ratio") + for item in csdm_object.x + if item.origin_offset != 0 + ] + data = csdm_object.y[0].components[0] + print(data) + # csdm_object.plot() + extra_dims = 1 + if len(csdm_object.x) > 2: + extra_dims = np.sum([item.coordinates.size for item in csdm_object.x[2:]]) + data.shape = (extra_dims, data.shape[-2], data.shape[-1]) + + reg_cq, reg_eta = (csdm_object.x[i].coordinates.value for i in range(2)) + d_cq = reg_cq[1] - reg_cq[0] + d_eta = reg_eta[1] - reg_eta[0] + sol = np.zeros((extra_dims, xq.count, yq.count)) + + bins = [xq.count, yq.count] + dx = xq.increment.value / 2 + dy = yq.increment.value / 2 + range_ = [ + [xq.coordinates[0].value - dx, xq.coordinates[-1].value + dx], + [yq.coordinates[0].value - dy, yq.coordinates[-1].value + dy], + ] + # print(range_) + avg_range_cq = (np.arange(n) - (n - 1) / 2) * d_cq / n + avg_range_eta = (np.arange(n) - (n - 1) / 2) * d_eta / n + # print(avg_range_cq) + # print(avg_range_eta) + for i, cq_item in enumerate(avg_range_cq): + for j, eta_item in enumerate(avg_range_eta): + # print(i,j) + cq__ = reg_cq + cq_item + # print(np.abs(cq__).shape) + eta__ = np.abs(reg_eta + eta_item) + cq_, eta_ = np.meshgrid(cq__, eta__) + # print(cq_) + # print() + cq_ = cq_.ravel() + eta_ = eta_.ravel() + # print(cq_.shape) + # print(eta_.shape) + # print(eta_) + # print() + + theta = np.ones(eta_.shape) * 1e10 + # print(theta.shape) + + index = np.where(cq_ > 0) + # print(f'pos index: {index}') + theta[index] = np.pi / 2 * (1 - eta_[index] / 2.0) + + index = np.where(cq_ <= 0) + # print(f'neg index: {index}') + + theta[index] = np.pi / 4.0 * eta_[index] + # print(cq_) + # print(np.where(theta * 180/np.pi==1e10)) + # print(list(zip(cq_, theta * 180/np.pi))) + + xq_grid = np.abs(cq_) * np.sin(theta) + yq_grid = np.abs(cq_) * np.cos(theta) + + # index = np.arange(xq.count) + # print(f'index: {index}') + # print(data[0].ravel()) + # print(data[0].shape) + ################eta_grid = (2 / np.pi) * np.arctan(y_ / x_) + # print(f'extra dims: {extra_dims}') + for i in range(extra_dims): + weight = deepcopy(data[i]).ravel() + # weight[index] /= 2 + # np.append(weight, weight[index]) + # plt.plot(weight) + # plt.show() + # print(np.where(weight>0.1)) + # print(range_) + sol_, _, _ = np.histogram2d( + xq_grid, yq_grid, weights=weight, bins=bins, range=range_ + ) + sol[i] += sol_ + # print(sol) + # print(np.where(sol>0.1)) + # plt.imshow(sol[0,:,:]) + # plt.imshow(sol_) + # print(xq.coordinates.shape) + # print(yq.coordinates.shape) + # plt.plot(weight) + # plt.show() + + sol /= n * n + # plt.plot(weight) + # plt.show() + del xq_grid, yq_grid, index, cq_, eta_, avg_range_cq, avg_range_eta + csdm_new = cp.as_csdm(np.squeeze(sol)) + csdm_new.x[0] = yq + csdm_new.x[1] = xq + if len(csdm_object.x) > 2: + csdm_new.x[2] = csdm_object.x[2] + return csdm_new + + +def get_polar_grids(ax, ticks=None, offset=0): + """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. + + Args: + ax: Matplotlib Axes. + ticks: Tick coordinates where radial grids are drawn. The value can be a list + or a numpy array. The default value is None. + offset: The grid is drawn at an offset away from the origin. + """ + ylim = ax.get_ylim() + xlim = ax.get_xlim() + if ticks is None: + x = np.asarray(ax.get_xticks()) + inc = x[1] - x[0] + size = x.size + x = np.arange(size + 5) * inc + else: + x = np.asarray(ticks) + + lw = 0.3 + t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) + t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) + + ax.add_artist(t1) + ax.add_artist(t2) + for x_ in x: + if x_ - offset > 0: + ax.add_artist( + plt.Circle( + (0, 0), + x_ - offset, + fill=False, + color="k", + linestyle="--", + linewidth=lw, + alpha=0.5, + ) + ) + + angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 4.0) + angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) + for ang_ in angle1: + ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) + for ang_ in angle2: + ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) + ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + +def get_quadpolar_grids(ax, ticks=None, offset=0): + """Generate a piece-wise polar grid of Haeberlen parameters, zeta and eta. + + Args: + ax: Matplotlib Axes. + ticks: Tick coordinates where radial grids are drawn. The value can be a list + or a numpy array. The default value is None. + offset: The grid is drawn at an offset away from the origin. + """ + ylim = ax.get_ylim() + xlim = ax.get_xlim() + if ticks is None: + x = np.asarray(ax.get_xticks()) + inc = x[1] - x[0] + size = x.size + x = np.arange(size + 5) * inc + else: + x = np.asarray(ticks) + + lw = 0.3 + # t1 = plt.Polygon([[0, 0], [0, x[-1]], [x[-1], x[-1]]], color="b", alpha=0.05) + # t2 = plt.Polygon([[0, 0], [x[-1], 0], [x[-1], x[-1]]], color="r", alpha=0.05) + + # ax.add_artist(t1) + # ax.add_artist(t2) + for x_ in x: + if x_ - offset > 0: + ax.add_artist( + plt.Circle( + (0, 0), + x_ - offset, + fill=False, + color="k", + linestyle="--", + linewidth=lw, + alpha=0.5, + ) + ) + + angle1 = np.tan(np.pi * np.asarray([0, 0.2, 0.4, 0.6, 0.8]) / 2.0) + # angle2 = np.tan(np.pi * np.asarray([0.8, 0.6, 0.4, 0.2, 0]) / 4.0) + for ang_ in angle1: + ax.plot(x, ((x - offset) * ang_) + offset, "k--", alpha=0.5, linewidth=lw) + # for ang_ in angle2: + # ax.plot(((x - offset) * ang_) + offset, x, "k--", alpha=0.5, linewidth=lw) + # ax.plot(x, x, "k", alpha=0.5, linewidth=2 * lw) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + +def plot_3d( + ax, + csdm_objects, + elev=28, + azim=-150, + x_lim=None, + y_lim=None, + z_lim=None, + cmap=cm.PiYG, + box=False, + clip_percent=0.0, + linewidth=1, + alpha=0.15, + **kwargs, +): + r"""Generate a 3D density plot with 2D contour and 1D projections. + + Args: + ax: Matplotlib Axes to render the plot. + csdm_objects: A 3D{1} CSDM object or a list of CSDM objects holding the data. + elev: (optional) The 3D view angle, elevation angle in the z plane. + azim: (optional) The 3D view angle, azimuth angle in the x-y plane. + x_lim: (optional) The x limit given as a list, [x_min, x_max]. + y_lim: (optional) The y limit given as a list, [y_min, y_max]. + z_lim: (optional) The z limit given as a list, [z_min, z_max]. + max_2d: (Optional) The normalization factor of the 2D contour projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the + maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. + max_1d: (Optional) The normalization factor of the 1D projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the + maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. + cmap: (Optional) The colormap or a list of colormaps used in rendering the + volumetric plot. The colormap list is applied to the ordered list of csdm_objects. + The same colormap is used for the 2D contour projections. For 1D plots, the first + color in the colormap scheme is used for the line color. + box: (Optional) If True, draw a box around the 3D data region. + clip_percent: (Optional) The amplitudes of the dataset below the given percent + is made transparent for the volumetric plot. + linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. + alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D + volume. + """ + csdm_object_list = csdm_objects + if not isinstance(csdm_objects, list): + csdm_object_list = [csdm_objects] + + if not isinstance(cmap, list): + cmap = [cmap] + + fraction = np.array([abs(item).sum() for item in csdm_object_list]) + fraction /= fraction.max() + + for index, item in enumerate(csdm_object_list): + plot_3d_( + ax, + item, + elev=elev, + azim=azim, + x_lim=x_lim, + y_lim=y_lim, + z_lim=z_lim, + fraction=fraction[index], + cmap=cmap[index], + box=box, + clip_percent=clip_percent, + linewidth=linewidth, + alpha=alpha, + **kwargs, + ) + + +def plot_3d_( + ax, + csdm_object, + elev=28, + azim=-150, + x_lim=None, + y_lim=None, + z_lim=None, + fraction=1.0, + cmap=cm.PiYG, + box=False, + clip_percent=0.0, + linewidth=1, + alpha=0.15, + **kwargs, +): + r"""Generate a 3D density plot with 2D contour and 1D projections. + + Args: + ax: Matplotlib Axes to render the plot. + csdm_object: A 3D{1} CSDM object holding the data. + elev: (optional) The 3D view angle, elevation angle in the z plane. + azim: (optional) The 3D view angle, azimuth angle in the x-y plane. + x_lim: (optional) The x limit given as a list, [x_min, x_max]. + y_lim: (optional) The y limit given as a list, [y_min, y_max]. + z_lim: (optional) The z limit given as a list, [z_min, z_max]. + max_2d: (Optional) The normalization factor of the 2D contour projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`yz`, `xz`, `xy`], where `ij` is the + maximum of the projection onto the `ij` plane, :math:`i,j \in [x, y, z]`. + max_1d: (Optional) The normalization factor of the 1D projections. The + attribute is meaningful when multiple 3D datasets are viewed on the same + plot. The value is given as a list, [`x`, `y`, `z`], where `i` is the + maximum of the projection onto the `i` axis, :math:`i \in [x, y, z]`. + cmap: (Optional) The colormap used in rendering the volumetric plot. The same + colormap is used for the 2D contour projections. For 1D plots, the first + color in the colormap scheme is used for the line color. + box: (Optional) If True, draw a box around the 3D data region. + clip_percent: (Optional) The amplitudes of the dataset below the given percent + is made transparent for the volumetric plot. + linewidth: (Optional) The linewidth of the 2D contours, 1D plots and box. + alpha: (Optional) The amount of alpha(transparency) applied in rendering the 3D + volume. + """ + lw = linewidth + if isinstance(csdm_object, cp.CSDM): + f = csdm_object.y[0].components[0].T + + a_, b_, c_ = (item for item in csdm_object.x) + + a = a_.coordinates.value + b = b_.coordinates.value + c = c_.coordinates.value + + xlabel = f"{a_.axis_label} - 0" + ylabel = f"{b_.axis_label} - 1" + zlabel = f"{c_.axis_label} - 2" + + else: + f = csdm_object + a = np.arange(f.shape[0]) + b = np.arange(f.shape[1]) + c = np.arange(f.shape[2]) + + xlabel = "x" + ylabel = "y" + zlabel = "z" + + delta_a = a[1] - a[0] + delta_b = b[1] - b[0] + delta_c = c[1] - c[0] + + clr = cmap + ck = cmap(0) + facecolors = cmap(f) + facecolors[:, :, :, -1] = f * alpha + index = np.where(f < clip_percent / 100) + facecolors[:, :, :, -1][index] = 0 + facecolors.shape = (f.size, 4) + + if x_lim is None: + x_lim = [a[0], a[-1]] + if y_lim is None: + y_lim = [b[0], b[-1]] + if z_lim is None: + z_lim = [c[0], c[-1]] + + offset_scalar = 50 + height_scalar = 10 + z_offset = (z_lim[1] - z_lim[0]) / offset_scalar + z_scale = np.abs(z_lim[1] - z_lim[0]) / height_scalar + offz_n = z_lim[0] - (delta_c / 2.0) + offz_1d = z_lim[1] + z_offset + + y_offset = (y_lim[1] - y_lim[0]) / offset_scalar + y_scale = np.abs(y_lim[1] - y_lim[0]) / height_scalar + offy = y_lim[1] + (delta_b / 2.0) + offy_n_1d = y_lim[0] - y_offset + + offx = x_lim[1] + (delta_a / 2.0) + x_scale = np.abs(x_lim[1] - x_lim[0]) / height_scalar + + if azim > -90 and azim <= 0: + offx = x_lim[0] - (delta_a / 2.0) + offy_n_1d = y_lim[0] - y_offset + + if azim > 0 and azim <= 90: + offy = y_lim[0] - (delta_b / 2.0) + offy_n_1d = y_lim[1] + y_offset + offx = x_lim[0] - (delta_a / 2.0) + + if azim > 90: + offx = x_lim[1] + (delta_a / 2.0) + offy_n_1d = y_lim[1] + y_offset + offy = y_lim[0] - (delta_b / 2.0) + + ax.set_proj_type("persp") + ax.view_init(elev=elev, azim=azim) + + # 2D x-y contour projection --------------------- + levels = (np.arange(20) + 1) / 20 + x1, y1 = np.meshgrid(a, b, indexing="ij") + dist_xy = f.mean(axis=2) + dist_xy *= fraction * z_scale / np.abs(dist_xy.max()) + ax.contour( + x1, + y1, + dist_xy, + zdir="z", + offset=offz_n, + cmap=clr, + levels=z_scale * levels, + linewidths=lw, + ) + + # 2D x-z contour projection + x1, y1 = np.meshgrid(a, c, indexing="ij") + dist_xz = f.mean(axis=1) + dist_xz *= fraction * y_scale / np.abs(dist_xz.max()) + ax.contour( + x1, + dist_xz, + y1, + zdir="y", + offset=offy, + cmap=clr, + levels=y_scale * levels, + linewidths=lw, + ) + + # 1D x-axis projection from 2D x-z projection + proj_x = dist_xz.mean(axis=1) + proj_x *= fraction * z_scale / np.abs(proj_x.max()) + ax.plot(a, np.sign(z_offset) * proj_x + offz_1d, offy, zdir="y", c=ck, linewidth=lw) + + # 1D z-axis projection from 2D x-z projection + proj_z = dist_xz.mean(axis=0) + proj_z *= fraction * y_scale / np.abs(proj_z.max()) + ax.plot( + np.sign(azim) * np.sign(y_offset) * proj_z + offy_n_1d, + c, + offx, + zdir="x", + c=ck, + linewidth=lw, + ) + ax.set_xlim(z_lim) + + # 2D y-z contour projection + x1, y1 = np.meshgrid(b, c, indexing="ij") + dist_yz = f.mean(axis=0) + dist_yz *= fraction * x_scale / np.abs(dist_yz.max()) + ax.contour( + dist_yz, + x1, + y1, + zdir="x", + offset=offx, + cmap=clr, + levels=x_scale * levels, + linewidths=lw, + ) + + # 1D y-axis projection + proj_y = dist_yz.mean(axis=1) + proj_y *= fraction * z_scale / np.abs(proj_y.max()) + ax.plot(b, np.sign(z_offset) * proj_y + offz_1d, offx, zdir="x", c=ck, linewidth=lw) + + ax.set_xlim(x_lim) + ax.set_ylim(y_lim) + ax.set_zlim(z_lim) + + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_zlabel(zlabel) + + x, y, z = np.meshgrid(a, b, c, indexing="ij") + + if "s" not in kwargs: + kwargs["s"] = 300 + ax.scatter(x.flat, y.flat, z.flat, marker="X", c=facecolors, **kwargs) + + # full box + r1 = [x_lim[0] - delta_a / 2, x_lim[-1] + delta_a / 2] + r2 = [y_lim[0] - delta_b / 2, y_lim[-1] + delta_b / 2] + r3 = [z_lim[-1] - delta_c / 2, z_lim[0] + delta_c / 2] + + l_box = lw + for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): + if np.sum(np.abs(s - e)) == r1[1] - r1[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + if np.sum(np.abs(s - e)) == r2[1] - r2[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + if np.sum(np.abs(s - e)) == r3[1] - r3[0]: + ax.plot3D(*zip(s, e), color="gray", linewidth=l_box) + + # draw cube + if box: + r1 = [a[0] - delta_a / 2, a[-1] + delta_a / 2] + r2 = [b[0] - delta_b / 2, b[-1] + delta_b / 2] + r3 = [c[0] - delta_c / 2, c[-1] + delta_c / 2] + + for s, e in combinations(np.array(list(product(r1, r2, r3))), 2): + if np.sum(np.abs(s - e)) == r1[1] - r1[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + if np.sum(np.abs(s - e)) == r2[1] - r2[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + if np.sum(np.abs(s - e)) == r3[1] - r3[0]: + ax.plot3D(*zip(s, e), c="blue", linestyle="dashed", linewidth=l_box) + + ax.xaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) + ax.yaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) + ax.zaxis._axinfo["grid"].update({"linewidth": 0.25, "color": "gray"}) diff --git a/setup.cfg b/setup.cfg index 0fbf5f6..a567d4a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,19 +1,19 @@ -[aliases] -test=pytest - -[tool:pytest] -addopts = - --ignore=docs/conf.py - --ignore=docs/galley_examples - --ignore=docs/_build - --ignore=examples - --ignore=setup.py - --doctest-glob='docs/*.rst' - --cov='./' - --doctest-modules - -[coverage:run] -omit = - setup.py - mrinversion/utils.py - mrinversion/linear_model/fista/__init__.py +[aliases] +test=pytest + +[tool:pytest] +addopts = + --ignore=docs/conf.py + --ignore=docs/galley_examples + --ignore=docs/_build + --ignore=examples + --ignore=setup.py + --doctest-glob='docs/*.rst' + --cov='./' + --doctest-modules + +[coverage:run] +omit = + setup.py + mrinversion/utils.py + mrinversion/linear_model/fista/__init__.py diff --git a/setup.py b/setup.py index c4432a3..004b563 100644 --- a/setup.py +++ b/setup.py @@ -1,68 +1,68 @@ -"""Setup for mrinversion package.""" -from os.path import abspath -from os.path import dirname -from os.path import join -from setuptools import find_packages, setup - - -with open("mrinversion/__init__.py", encoding="utf-8") as f: - for line in f.readlines(): - if "__version__" in line: - before_keyword, keyword, after_keyword = line.partition("=") - version = after_keyword.strip()[1:-1] - -module_dir = dirname(abspath(__file__)) - -install_requires = [ - "numpy>=2.0", - "setuptools>=27.3", - "csdmpy>=0.7", - "mrsimulator>=1.0.0", - "scikit-learn>=1.5.2", - "numba>=0.61.2", -] - -setup_requires = ["setuptools>=27.3", "numpy"] -extras = {"matplotlib": ["matplotlib>=3.0"]} - - -setup( - name="mrinversion", - version=version, - description=( - "Python based statistical learning of NMR tensor and relaxation parameters " - "distribution." - ), - long_description=open(join(module_dir, "README.md"), encoding="utf-8").read(), - long_description_content_type="text/markdown", - author="Deepansh J. Srivastava", - author_email="deepansh2012@gmail.com", - python_requires=">=3.10", - url="https://github.com/DeepanshS/mrinversion/", - packages=find_packages(), - install_requires=install_requires, - setup_requires=setup_requires, - extras_require=extras, - tests_require=["pytest", "pytest-runner"], - include_package_data=True, - zip_safe=False, - license="BSD-3-Clause", - classifiers=[ - # Trove classifiers - # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers - "Intended Audience :: Science/Research", - "Intended Audience :: Developers", - "Operating System :: OS Independent", - "Development Status :: 4 - Beta", - "License :: OSI Approved :: BSD License", - "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.10", - "Programming Language :: Python :: 3.11", - "Programming Language :: Python :: 3.12", - "Programming Language :: Python :: 3.13", - "Topic :: Scientific/Engineering", - "Topic :: Education", - "Topic :: Scientific/Engineering :: Chemistry", - "Topic :: Scientific/Engineering :: Physics", - ], -) +"""Setup for mrinversion package.""" +from os.path import abspath +from os.path import dirname +from os.path import join +from setuptools import find_packages, setup + + +with open("mrinversion/__init__.py", encoding="utf-8") as f: + for line in f.readlines(): + if "__version__" in line: + before_keyword, keyword, after_keyword = line.partition("=") + version = after_keyword.strip()[1:-1] + +module_dir = dirname(abspath(__file__)) + +install_requires = [ + "numpy>=2.0", + "setuptools>=27.3", + "csdmpy>=0.7", + "mrsimulator>=1.0.0", + "scikit-learn>=1.5.2", + "numba>=0.61.2", +] + +setup_requires = ["setuptools>=27.3", "numpy"] +extras = {"matplotlib": ["matplotlib>=3.0"]} + + +setup( + name="mrinversion", + version=version, + description=( + "Python based statistical learning of NMR tensor and relaxation parameters " + "distribution." + ), + long_description=open(join(module_dir, "README.md"), encoding="utf-8").read(), + long_description_content_type="text/markdown", + author="Deepansh J. Srivastava", + author_email="deepansh2012@gmail.com", + python_requires=">=3.10", + url="https://github.com/DeepanshS/mrinversion/", + packages=find_packages(), + install_requires=install_requires, + setup_requires=setup_requires, + extras_require=extras, + tests_require=["pytest", "pytest-runner"], + include_package_data=True, + zip_safe=False, + license="BSD-3-Clause", + classifiers=[ + # Trove classifiers + # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers + "Intended Audience :: Science/Research", + "Intended Audience :: Developers", + "Operating System :: OS Independent", + "Development Status :: 4 - Beta", + "License :: OSI Approved :: BSD License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Topic :: Scientific/Engineering", + "Topic :: Education", + "Topic :: Scientific/Engineering :: Chemistry", + "Topic :: Scientific/Engineering :: Physics", + ], +)