From 73096872253f6e5e42dcf62c27f838a841246e77 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 06:57:15 -0400 Subject: [PATCH 01/27] numpy<2.0 --- environment-dev.yml | 2 +- environment-docs.yml | 2 +- environment.yml | 2 +- requirements.txt | 2 +- setup.py | 6 +++--- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/environment-dev.yml b/environment-dev.yml index 71acebe27..b0193b7d1 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -5,7 +5,7 @@ channels: dependencies: - fftw>=3.3.9 - openblas - - numpy>=1.20 + - numpy<2.0 - nomkl - setuptools>=27.3 - cython>=0.29 diff --git a/environment-docs.yml b/environment-docs.yml index 4a6207581..798e1fb4e 100644 --- a/environment-docs.yml +++ b/environment-docs.yml @@ -5,7 +5,7 @@ channels: dependencies: - fftw>=3.3.9 - openblas - - numpy>=1.20 + - numpy<2.0 - nomkl - setuptools>=27.3 - cython>=0.29 diff --git a/environment.yml b/environment.yml index 101588c49..df559f223 100644 --- a/environment.yml +++ b/environment.yml @@ -10,7 +10,7 @@ dependencies: - setuptools>=27.3 - pip - pip: - - numpy>=1.20 + - numpy<2.0 - cython>=0.29 - matplotlib>=3.3.4 - joblib>=1.0.0 diff --git a/requirements.txt b/requirements.txt index 311bc6f57..937a9e415 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ -numpy>=1.20 +numpy<2.0 matplotlib>=3.3.4 csdmpy>=0.6 astropy>=5.1 diff --git a/setup.py b/setup.py index 496f55d5b..7e82cfd1b 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ def message(lib, env, command, key): "For example,\n", '\texport LDFLAGS="-L/usr/local/opt/openblas/lib"\n', '\texport CPPFLAGS="-I/usr/local/opt/openblas/include"\n', - f"\nYou can also try installing '{lib}' from {env} with:", + f"\nYou can also try installing '{lib}' from {env} with: ", f"\n\t{command} install {arg}\n", ) warnings.warn("".join(warning)) @@ -368,9 +368,9 @@ def fftw_info(self): url="https://github.com/deepanshs/mrsimulator/", packages=find_packages("src"), package_dir={"": "src"}, - setup_requires=["numpy>=1.20"], + setup_requires=["numpy<2.0"], install_requires=[ - "numpy>=1.20", + "numpy<2.0", "csdmpy>=0.6", "pydantic>=2.8", "monty>=2.0.4", From 1ed871aa914215b3ae83b6e08479d3cce1864fee Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 06:57:35 -0400 Subject: [PATCH 02/27] numpy<2.0_Re --- docs/installation/requirements.rst | 2 +- docs/requirements.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/installation/requirements.rst b/docs/installation/requirements.rst index 9ef64dcc7..8088213d2 100644 --- a/docs/installation/requirements.rst +++ b/docs/installation/requirements.rst @@ -16,7 +16,7 @@ following operating systems: **Required packages** -- `numpy>=1.20 `_ +- `numpy<2.0 `_ - `matplotlib>=3.3.4 `_ for figures and visualization - `lmfit>=1.0.3 `_ for least-squares fitting - `pandas>=1.1.3 `_ diff --git a/docs/requirements.txt b/docs/requirements.txt index 7e7a1d01d..ffa875a77 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,4 +1,4 @@ -numpy>=1.20 +numpy<2.0 matplotlib>=3.3.3 csdmpy>=0.6 astropy>=5.1 From 09066592d857b6d9646e33260abb507ebc5f4a14 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 21:55:11 -0400 Subject: [PATCH 03/27] force gcc compiler for mac --- .github/workflows/continuous-integration-pip.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index cb512b19b..5ad743236 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -100,7 +100,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: python setup.py develop + run: CC=gcc python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -134,7 +134,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: python setup.py develop + run: CC=gcc python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml From 0671369fba91f74e32a2b3c77c10bcfbe2b769d7 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:08:02 -0400 Subject: [PATCH 04/27] force gcc-14 compiler for mac --- .github/workflows/continuous-integration-pip.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 5ad743236..a08c8c935 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -100,7 +100,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc python setup.py develop + run: CC=gcc-14 python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -134,7 +134,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc python setup.py develop + run: CC=gcc-14 python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml From 7a40f9a1d6cc900b93ad5c1d37b189766e3143ba Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:14:22 -0400 Subject: [PATCH 05/27] update mac setup --- setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 7e82cfd1b..589e89b36 100644 --- a/setup.py +++ b/setup.py @@ -183,8 +183,9 @@ def __init__(self): # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", # "-Rpass-analysis=loop-vectorize", - "-fvectorize", + # "-fvectorize", "-fcommon", + "-Wno-incompatible-pointer-types", ] self.extra_link_args += ["-lm"] From ead7beb0280ec6e1aac0e24f57dd8559c65d398b Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:19:49 -0400 Subject: [PATCH 06/27] O2 compile --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 589e89b36..17413dad6 100644 --- a/setup.py +++ b/setup.py @@ -178,7 +178,7 @@ class MacOSSetup(Setup): def __init__(self): super().__init__() self.extra_compile_args = [ - "-O3", + # "-O3", "-ffast-math", # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", From c5b25643c41df88fba729934f946248eccbc09ac Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:22:49 -0400 Subject: [PATCH 07/27] O2 compile --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 17413dad6..31977fa07 100644 --- a/setup.py +++ b/setup.py @@ -178,7 +178,7 @@ class MacOSSetup(Setup): def __init__(self): super().__init__() self.extra_compile_args = [ - # "-O3", + "-O2", "-ffast-math", # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", From a9cffc0690699302ce06473421a75c43fb5d4a05 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:34:58 -0400 Subject: [PATCH 08/27] remove bit operation --- src/c_lib/include/vm.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/c_lib/include/vm.h b/src/c_lib/include/vm.h index f1719174e..020c477d1 100644 --- a/src/c_lib/include/vm.h +++ b/src/c_lib/include/vm.h @@ -18,8 +18,9 @@ // absolute double value static inline double absd(double a) { - *((unsigned __int64_ *)&a) &= ~(1ULL << 63); - return a; + // *((unsigned __int64_ *)&a) &= ~(1ULL << 63); + // return a; + return fabs(a); } /** Arithmetic suit ======================================================== */ From 99f0562398fe5f8f7107aaa502259e65095dce85 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:41:14 -0400 Subject: [PATCH 09/27] . --- .github/workflows/continuous-integration-pip.yml | 4 ++-- setup.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index a08c8c935..cb512b19b 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -100,7 +100,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc-14 python setup.py develop + run: python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -134,7 +134,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc-14 python setup.py develop + run: python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml diff --git a/setup.py b/setup.py index 31977fa07..a1b473eff 100644 --- a/setup.py +++ b/setup.py @@ -178,14 +178,14 @@ class MacOSSetup(Setup): def __init__(self): super().__init__() self.extra_compile_args = [ - "-O2", + "-O1", "-ffast-math", # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", # "-Rpass-analysis=loop-vectorize", # "-fvectorize", "-fcommon", - "-Wno-incompatible-pointer-types", + # "-Wno-incompatible-pointer-types", ] self.extra_link_args += ["-lm"] From 928210c4caae6389a7d649d343be5159a972641c Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Mon, 19 Aug 2024 22:46:41 -0400 Subject: [PATCH 10/27] . --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index a1b473eff..19980f4c0 100644 --- a/setup.py +++ b/setup.py @@ -179,7 +179,7 @@ def __init__(self): super().__init__() self.extra_compile_args = [ "-O1", - "-ffast-math", + # "-ffast-math", # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", # "-Rpass-analysis=loop-vectorize", From 7d9ad6ef137e88f97c114aec7868cf5789bcd4aa Mon Sep 17 00:00:00 2001 From: Deepansh Srivastava <21365911+deepanshs@users.noreply.github.com> Date: Tue, 20 Aug 2024 08:39:00 -0400 Subject: [PATCH 11/27] test --- .github/workflows/continuous-integration-pip.yml | 4 ++-- setup.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index cb512b19b..d1d07fa30 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -100,7 +100,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: python setup.py develop + run: CC=gcc-13 python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -134,7 +134,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: python setup.py develop + run: CC=gcc-12 python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml diff --git a/setup.py b/setup.py index 19980f4c0..282471bf3 100644 --- a/setup.py +++ b/setup.py @@ -185,7 +185,7 @@ def __init__(self): # "-Rpass-analysis=loop-vectorize", # "-fvectorize", "-fcommon", - # "-Wno-incompatible-pointer-types", + "-Wno-incompatible-pointer-types", ] self.extra_link_args += ["-lm"] From f063a7052b77e637c1b3addc340439c9ff651fd1 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 20 Aug 2024 17:36:28 -0400 Subject: [PATCH 12/27] . --- .../workflows/continuous-integration-pip.yml | 155 +++++++++--------- 1 file changed, 80 insertions(+), 75 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index d1d07fa30..05925a0f2 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -38,48 +38,49 @@ jobs: # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide flake8 . --count --exit-zero --max-complexity=8 --max-line-length=88 --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd" - testing_unix: - needs: [code_lint] - runs-on: "ubuntu-latest" - strategy: - matrix: - python-version: [3.9, "3.10", "3.11", "3.12"] - - steps: - - name: Checkout - uses: actions/checkout@v4 - - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v5 - with: - python-version: ${{ matrix.python-version }} - - - name: Install linux system dependencies - run: | - sudo apt-get install --yes libopenblas-dev libfftw3-dev - gcc -v - - - 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@v4.5.0 + # testing_unix: + # needs: [code_lint] + # runs-on: "ubuntu-latest" + # strategy: + # matrix: + # python-version: [3.9, "3.10", "3.11", "3.12"] + + # steps: + # - name: Checkout + # uses: actions/checkout@v4 + + # - name: Set up Python ${{ matrix.python-version }} + # uses: actions/setup-python@v5 + # with: + # python-version: ${{ matrix.python-version }} + + # - name: Install linux system dependencies + # run: | + # sudo apt-get install --yes libopenblas-dev libfftw3-dev + # gcc -v + + # - 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@v4.5.0 testing_mac_arm: needs: [code_lint] runs-on: "macos-latest" strategy: matrix: - python-version: [3.9, "3.10", "3.11", "3.12"] + python-version: ["3.10"] + # python-version: [3.9, "3.10", "3.11", "3.12"] steps: - name: Checkout @@ -91,7 +92,9 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install macos system dependencies - run: brew install openblas fftw + run: | + brew install openblas fftw + which clang - name: Install dependencies run: | @@ -100,7 +103,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc-13 python setup.py develop + run: python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -125,7 +128,9 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install macos system dependencies - run: brew install openblas fftw + run: | + brew install openblas fftw + which clang - name: Install dependencies run: | @@ -134,7 +139,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc-12 python setup.py develop + run: python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -142,38 +147,38 @@ jobs: - name: Upload coverage uses: codecov/codecov-action@v4.5.0 - testing_windows: - needs: [code_lint] - runs-on: "windows-latest" - strategy: - matrix: - python-version: [3.9, "3.10", "3.11", "3.12"] - - steps: - - name: Checkout - uses: actions/checkout@v4 - - name: Setup Miniconda - uses: conda-incubator/setup-miniconda@v3 - env: - ACTIONS_ALLOW_UNSECURE_COMMANDS: "true" - with: - auto-update-conda: true - auto-activate-base: false - miniconda-version: "latest" - python-version: ${{ matrix.python-version }} - environment-file: environment-dev.yml - activate-environment: mrsimulator-dev - - run: | - conda --version - which python - - name: Build and install package from source - shell: pwsh - run: | - conda --version - which python - python setup.py develop - - name: Test with pytest - shell: pwsh - run: pytest --cov=./ --cov-report=xml - - name: Upload coverage - uses: codecov/codecov-action@v4.5.0 + # testing_windows: + # needs: [code_lint] + # runs-on: "windows-latest" + # strategy: + # matrix: + # python-version: [3.9, "3.10", "3.11", "3.12"] + + # steps: + # - name: Checkout + # uses: actions/checkout@v4 + # - name: Setup Miniconda + # uses: conda-incubator/setup-miniconda@v3 + # env: + # ACTIONS_ALLOW_UNSECURE_COMMANDS: "true" + # with: + # auto-update-conda: true + # auto-activate-base: false + # miniconda-version: "latest" + # python-version: ${{ matrix.python-version }} + # environment-file: environment-dev.yml + # activate-environment: mrsimulator-dev + # - run: | + # conda --version + # which python + # - name: Build and install package from source + # shell: pwsh + # run: | + # conda --version + # which python + # python setup.py develop + # - name: Test with pytest + # shell: pwsh + # run: pytest --cov=./ --cov-report=xml + # - name: Upload coverage + # uses: codecov/codecov-action@v4.5.0 From 019d87941c57cb834d937045e8b8cb88ff7deb40 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 20 Aug 2024 17:38:59 -0400 Subject: [PATCH 13/27] . --- .github/workflows/continuous-integration-pip.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 05925a0f2..509287681 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -95,6 +95,7 @@ jobs: run: | brew install openblas fftw which clang + clang --version - name: Install dependencies run: | @@ -131,6 +132,7 @@ jobs: run: | brew install openblas fftw which clang + clang --version - name: Install dependencies run: | From a2615716722c70ad921346a826bac4da8d30bd70 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 20 Aug 2024 17:53:58 -0400 Subject: [PATCH 14/27] . --- .github/workflows/continuous-integration-pip.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 509287681..9b4e324cc 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -141,7 +141,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: python setup.py develop + run: CC=gcc-12 python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml From be1e3a7b22ad28daee2dc0588daa7a12089fa543 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 20 Aug 2024 18:13:55 -0400 Subject: [PATCH 15/27] . --- docs/installation/requirements.rst | 2 +- docs/requirements.txt | 2 +- environment-dev.yml | 2 +- environment-docs.yml | 2 +- environment.yml | 2 +- requirements.txt | 2 +- setup.py | 9 ++++----- 7 files changed, 10 insertions(+), 11 deletions(-) diff --git a/docs/installation/requirements.rst b/docs/installation/requirements.rst index 8088213d2..9ef64dcc7 100644 --- a/docs/installation/requirements.rst +++ b/docs/installation/requirements.rst @@ -16,7 +16,7 @@ following operating systems: **Required packages** -- `numpy<2.0 `_ +- `numpy>=1.20 `_ - `matplotlib>=3.3.4 `_ for figures and visualization - `lmfit>=1.0.3 `_ for least-squares fitting - `pandas>=1.1.3 `_ diff --git a/docs/requirements.txt b/docs/requirements.txt index ffa875a77..7e7a1d01d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,4 +1,4 @@ -numpy<2.0 +numpy>=1.20 matplotlib>=3.3.3 csdmpy>=0.6 astropy>=5.1 diff --git a/environment-dev.yml b/environment-dev.yml index b0193b7d1..71acebe27 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -5,7 +5,7 @@ channels: dependencies: - fftw>=3.3.9 - openblas - - numpy<2.0 + - numpy>=1.20 - nomkl - setuptools>=27.3 - cython>=0.29 diff --git a/environment-docs.yml b/environment-docs.yml index 798e1fb4e..4a6207581 100644 --- a/environment-docs.yml +++ b/environment-docs.yml @@ -5,7 +5,7 @@ channels: dependencies: - fftw>=3.3.9 - openblas - - numpy<2.0 + - numpy>=1.20 - nomkl - setuptools>=27.3 - cython>=0.29 diff --git a/environment.yml b/environment.yml index df559f223..101588c49 100644 --- a/environment.yml +++ b/environment.yml @@ -10,7 +10,7 @@ dependencies: - setuptools>=27.3 - pip - pip: - - numpy<2.0 + - numpy>=1.20 - cython>=0.29 - matplotlib>=3.3.4 - joblib>=1.0.0 diff --git a/requirements.txt b/requirements.txt index 937a9e415..311bc6f57 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ -numpy<2.0 +numpy>=1.20 matplotlib>=3.3.4 csdmpy>=0.6 astropy>=5.1 diff --git a/setup.py b/setup.py index 282471bf3..2e8c8d363 100644 --- a/setup.py +++ b/setup.py @@ -178,14 +178,13 @@ class MacOSSetup(Setup): def __init__(self): super().__init__() self.extra_compile_args = [ - "-O1", - # "-ffast-math", + "-O3", + "-ffast-math", # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", # "-Rpass-analysis=loop-vectorize", # "-fvectorize", "-fcommon", - "-Wno-incompatible-pointer-types", ] self.extra_link_args += ["-lm"] @@ -369,9 +368,9 @@ def fftw_info(self): url="https://github.com/deepanshs/mrsimulator/", packages=find_packages("src"), package_dir={"": "src"}, - setup_requires=["numpy<2.0"], + setup_requires=["numpy>=1.20"], install_requires=[ - "numpy<2.0", + "numpy>=1.20", "csdmpy>=0.6", "pydantic>=2.8", "monty>=2.0.4", From ef95d9769ac75797e95f1f16b8bbc839574aab8e Mon Sep 17 00:00:00 2001 From: unknown <21365911+deepanshs@users.noreply.github.com> Date: Sun, 25 Aug 2024 15:41:39 -0400 Subject: [PATCH 16/27] apply compile warning suggestions --- src/c_lib/base/base_model.pyx | 36 +++++++++++++------------- src/c_lib/clib/clib.pyx | 4 +-- src/c_lib/include/config.h | 5 ++-- src/c_lib/lib/interpolation.c | 37 ++++++++++++++------------- src/c_lib/lib/method.c | 32 ++++++++++++----------- src/c_lib/test/test.pyx | 48 +++++++++++++++++------------------ 6 files changed, 84 insertions(+), 78 deletions(-) diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index 612f72f4c..f2edb9140 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -64,7 +64,7 @@ def core_simulator(method, else: positions = None averaging_scheme = clib.MRS_create_averaging_scheme_from_alpha_beta( - alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=alpha.size, + alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=np.uint32(alpha.size), allow_4th_rank=allow_4th_rank, n_gamma=number_of_gamma_angles, position_size=position_size, positions=&positions[0], interpolation=interpolation ) @@ -76,7 +76,7 @@ def core_simulator(method, ) # create C spectral dimensions ________________________________________________ - cdef int n_dimension = len(method.spectral_dimensions) + cdef int n_dimension = int(len(method.spectral_dimensions)) # Define and allocate C numpy arrays n_points = 1 @@ -257,7 +257,7 @@ def core_simulator(method, # sub_sites = [site for site in spin_sys.sites if site.isotope.symbol == isotope] # index_.append(index) - number_of_sites = len(spin_sys.sites) + number_of_sites = int(len(spin_sys.sites)) # ------------------------------------------------------------------------ # Site specification @@ -324,7 +324,7 @@ def core_simulator(method, if debug: print(f'Quadrupolar coupling constant (Cq) = {Cq_e[i]/1e6} MHz') print(f'Quadrupolar asymmetry (η) = {eta_e}') - print(f'Quadrupolar orientation (alpha/beta/gamma) = {ori_e}]') + print(f'Quadrupolar orientation (alpha/beta/gamma) = {ori_e}') # sites packed as c struct sites_c.number_of_sites = number_of_sites @@ -345,7 +345,7 @@ def core_simulator(method, # J-coupling couplings_c.number_of_couplings = 0 if spin_sys.couplings is not None: - number_of_couplings = len(spin_sys.couplings) + number_of_couplings = int(len(spin_sys.couplings)) spin_index_ij = np.zeros(2*number_of_couplings, dtype=np.int32) iso_j = np.zeros(number_of_couplings, dtype=np.float64) @@ -395,17 +395,17 @@ def core_simulator(method, if dipolar.gamma is not None: ori_d[i3+2] = dipolar.gamma - # if debug: - # print(f'N couplings = {number_of_couplings}') - # print(f'site index J = {spin_index_ij}') - # print(f'Isotropic J = {iso_j} Hz') - # print(f'J anisotropy = {zeta_j} Hz') - # print(f'J asymmetry = {eta_j}') - # print(f'J orientation = {ori_j}') + if debug: + print(f'N couplings = {number_of_couplings}') + print(f'site index J = {spin_index_ij}') + print(f'Isotropic J = {iso_j} Hz') + print(f'J anisotropy = {zeta_j} Hz') + print(f'J asymmetry = {eta_j}') + print(f'J orientation = {ori_j}') - # print(f'Dipolar coupling constant = {D_d} Hz') - # print(f'Dipolar asymmetry = {eta_d}') - # print(f'Dipolar orientation = {ori_d}') + print(f'Dipolar coupling constant = {D_d} Hz') + print(f'Dipolar asymmetry = {eta_d}') + print(f'Dipolar orientation = {ori_d}') # couplings packed as c struct couplings_c.number_of_couplings = number_of_couplings @@ -493,7 +493,7 @@ def core_simulator(method, @cython.boundscheck(False) @cython.wraparound(False) def get_zeeman_states(sys): - cdef int i, j, n_site = len(sys.sites) + cdef int i, j, n_site = int(len(sys.sites)) two_Ip1 = [int(2 * site.isotope.spin + 1) for site in sys.sites] spin_quantum_numbers = [ @@ -593,7 +593,7 @@ def calculate_transition_connect_weight( Return: A complex amplitude. """ cdef ndarray[double] factor = np.asarray([1, 0], dtype=np.float64) - cdef int i, n_sites = spin.size + cdef int i, n_sites = int(spin.size) cdef float m1_f, m1_i, m2_f, m2_i for i in range(n_sites): m1_f = trans1[1][i] # starting transition final state @@ -617,7 +617,7 @@ def get_Haeberlen_components( double eta, double rho): """Return random extended czjzek tensors in Haeberlen convention""" - cdef int n = expr_base_p.shape[1] + cdef int n = int(expr_base_p.shape[1]) cdef ndarray[double, ndim=2] param = np.empty((n, 2), dtype=float) clib.vm_haeberlen_components( n, &expr_base_p[0, 0], &expr_base_q[0, 0], zeta, eta, rho, ¶m[0, 0] diff --git a/src/c_lib/clib/clib.pyx b/src/c_lib/clib/clib.pyx index e42ef705f..8aeec6842 100644 --- a/src/c_lib/clib/clib.pyx +++ b/src/c_lib/clib/clib.pyx @@ -20,7 +20,7 @@ def histogram1d( ndarray[double, ndim=1] weights, ): cdef ndarray[double] hist = np.zeros(x_count, dtype=np.float64) - cdef int sample_count = sample_x.size + cdef int sample_count = int(sample_x.size) cdef int c_int64 = ctypes.sizeof(ctypes.c_int64) cdef int stride_x = sample_x.ctypes.strides[0] / c_int64 cdef int stride_w = weights.ctypes.strides[0] / c_int64 @@ -49,7 +49,7 @@ def histogram2d( int interp=0, ): cdef ndarray[double, ndim=2] hist = np.zeros((x_count, y_count), dtype=np.float64) - cdef int sample_count = sample_x.size + cdef int sample_count = int(sample_x.size) cdef int c_int64 = ctypes.sizeof(ctypes.c_int64) cdef int stride_x = sample_x.ctypes.strides[0] / c_int64 cdef int stride_y = sample_y.ctypes.strides[0] / c_int64 diff --git a/src/c_lib/include/config.h b/src/c_lib/include/config.h index cbdf45a2c..18f9ab54d 100644 --- a/src/c_lib/include/config.h +++ b/src/c_lib/include/config.h @@ -60,12 +60,13 @@ typedef float complex64[2]; // ---------------------------------------------------------------------------- // // OS base definitions -------------------------------------------------------- // -#define __int64_ long // for both mac-os and linux (unix system) - #ifdef _WIN32 // windows 32-bit or 64-bit #define __int64_ long long // #if _MSC_VER && !__INTEL_COMPILER // windows msvc compiler // #endif // end windows 32-bit or 64-bit +#else +// for both mac-os and linux (unix system) +#define __int64_ long #endif // end windows msvc compiler // ---------------------------------------------------------------------------- // diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index 633fd078a..15c31a984 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -72,10 +72,10 @@ static void inline delta_fn_gauss_interpolation(const double *freq, const int *p index = (int)res; w = res - (double)index; - ip2 = 2 * gauss_table_precision_inverse - index; - ip1 = gauss_table_precision_inverse - index; - im1 = gauss_table_precision_inverse + index; - im2 = 2 * gauss_table_precision_inverse + index; + ip2 = (int)(2 * gauss_table_precision_inverse - index); + ip1 = (int)(gauss_table_precision_inverse - index); + im1 = (int)(gauss_table_precision_inverse + index); + im2 = (int)(2 * gauss_table_precision_inverse + index); a0 = lerp_plus(w, im2); a1 = lerp_plus(w, im1); @@ -263,26 +263,29 @@ static inline void __triangle_interpolation(double *freq1, double *freq2, double void triangle_interpolation1D(double *freq1, double *freq2, double *freq3, double *amp, double *spec, int *points, unsigned int iso_intrp) { if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { - if (iso_intrp == 0) return delta_fn_linear_interpolation(freq1, points, amp, spec); - if (iso_intrp == 1) return delta_fn_gauss_interpolation(freq1, points, amp, spec); + if (iso_intrp == 0) delta_fn_linear_interpolation(freq1, points, amp, spec); + if (iso_intrp == 1) delta_fn_gauss_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); } - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); } void triangle_interpolation1D_linear(double *freq1, double *freq2, double *freq3, double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) - return delta_fn_linear_interpolation(freq1, points, amp, spec); - - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { + delta_fn_linear_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + } } void triangle_interpolation1D_gaussian(double *freq1, double *freq2, double *freq3, double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) - return delta_fn_gauss_interpolation(freq1, points, amp, spec); - - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { + delta_fn_gauss_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + } } /** @@ -838,8 +841,8 @@ void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *a int_i_stride += stride; int_j_stride += stride; } - if (iso_intrp == 0) return delta_fn_linear_interpolation(freq, &n_spec, &1, spec); - if (iso_intrp == 1) return delta_fn_gauss_interpolation(freq, &n_spec, &1, spec); + if (iso_intrp == 0) delta_fn_linear_interpolation(freq, &n_spec, &1, spec); + if (iso_intrp == 1) delta_fn_gauss_interpolation(freq, &n_spec, &1, spec); } void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, diff --git a/src/c_lib/lib/method.c b/src/c_lib/lib/method.c index 341b1c784..f521bf3d4 100644 --- a/src/c_lib/lib/method.c +++ b/src/c_lib/lib/method.c @@ -213,21 +213,23 @@ MRS_dimension *MRS_create_dimensions( rotor_angle_in_rad += n_events[i]; magnetic_flux_density_in_T += n_events[i]; - // printf("dimension %d\n", i); - // printf("\tcount %d\n", dimension[i].count); - // printf("\tincrement %f Hz\n", dimension[i].increment); - // printf("\tcoordinates offset %f Hz\n", dimension[i].coordinates_offset); - // for (int j = 0; j < n_events[i]; j++) { - // printf("\tEvent %d\n", j); - // printf("\t\tfraction %f\n", dimension[i].events[j].fraction); - // printf("\t\trotor frequency %f Hz\n", - // dimension[i].events[j].rotor_frequency_in_Hz); - // printf("\t\trotor angle %f rad\n", - // dimension[i].events[j].rotor_angle_in_rad); - // printf("\t\tmagnetic flux density %f T\n", - // dimension[i].events[j].magnetic_flux_density_in_T); - // printf("\n"); - // } + if (DEBUG) { + printf("Creating dimensions"); + printf("dimension %d\n", i); + printf("\tcount %d\n", dimension[i].count); + printf("\tincrement %f Hz\n", dimension[i].increment); + printf("\tcoordinates offset %f Hz\n", dimension[i].coordinates_offset); + for (int j = 0; j < n_events[i]; j++) { + printf("\tEvent %d\n", j); + printf("\t\tfraction %f\n", dimension[i].events[j].fraction); + printf("\t\trotor frequency %f Hz\n", + dimension[i].events[j].rotor_frequency_in_Hz); + printf("\t\trotor angle %f rad\n", dimension[i].events[j].rotor_angle_in_rad); + printf("\t\tmagnetic flux density %f T\n", + dimension[i].events[j].magnetic_flux_density_in_T); + printf("\n"); + } + } } return dimension; } diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index ba6bc1f11..7b076559c 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -23,7 +23,7 @@ def wigner_d_element(float l, float m1, float m2, double beta): @cython.boundscheck(False) @cython.wraparound(False) def wigner_d_matrices(int l, np.ndarray[double] angle): - cdef int n = angle.size + cdef int n = int(angle.size) cdef int n1 = (2*l+1)**2 cdef np.ndarray[double] wigner = np.empty(n*n1, dtype=np.float64) clib.wigner_d_matrices(l, n, &angle[0], &wigner[0]) @@ -46,7 +46,7 @@ def wigner_d_matrices_from_exp_I_beta(int l, bool_t half, np.ndarray[double comp :ivar exp_I_beta: An 1D numpy array or a scalar representing $\exp\beta$. """ n1 = (2 * l + 1) - cdef int n = exp_I_beta.size + cdef int n = int(exp_I_beta.size) size = n * n1*(l+1) if half else n * n1**2 cdef np.ndarray[double, ndim=1] wigner = np.empty(size) @@ -88,7 +88,7 @@ def __wigner_rotation_2(int l, np.ndarray[double] cos_alpha, np.ndarray[double complex] R_in): cdef int n1 = 2 * l + 1 - cdef int n = cos_alpha.size + cdef int n = int(cos_alpha.size) cdef np.ndarray[double, ndim=1] wigner cdef np.ndarray[double complex, ndim=1] exp_I_beta wigner = np.empty(n1 * (l+1) * n, dtype=np.float64) @@ -173,9 +173,9 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): @cython.wraparound(False) def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, int stride=1): cdef int i - cdef int number_of_sidebands = amp.shape[0] + cdef int number_of_sidebands = int(amp.shape[0]) for i in range(number_of_sidebands): - clib.octahedronInterpolation(&spec[0], &freq[i,0], nt, &[i,0], stride, spec.size) + clib.octahedronInterpolation(&spec[0], &freq[i,0], nt, &[i,0], stride, int(spec.size)) @cython.boundscheck(False) @@ -283,122 +283,122 @@ def vm_absd(double a): @cython.wraparound(False) def vm_add(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_add(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_add(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_add_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_add_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_add_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_sub(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_subtract(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_subtract(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sub_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_subtract_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_subtract_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_mult(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_multiply(A.size, &A[0], 1, &B[0], &res[0]) + clib.test_vm_double_multiply(int(A.size), &A[0], 1, &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_mult_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_multiply_inplace(A.size, &A[0], 1, &B[0], 1) + clib.test_vm_double_multiply_inplace(int(A.size), &A[0], 1, &B[0], 1) @cython.boundscheck(False) @cython.wraparound(False) def vm_div(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_divide(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_divide(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_div_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_divide_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_divide_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_cmult(np.ndarray[complex] A, np.ndarray[complex] B): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_multiply(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_complex_multiply(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_cmult_conj(np.ndarray[complex] A, np.ndarray[complex] B): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_conj_multiply(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_complex_conj_multiply(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sq(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square(A.size, &A[0], &res[0]) + clib.test_vm_double_square(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sq_inplace(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_inplace(A.size, &A[0]) + clib.test_vm_double_square_inplace(int(A.size), &A[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_sqrt(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_root(A.size, &A[0], &res[0]) + clib.test_vm_double_square_root(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sqrt_inplace(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_root_inplace(A.size, &A[0]) + clib.test_vm_double_square_root_inplace(int(A.size), &A[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_cos(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_cosine(A.size, &A[0], &res[0]) + clib.test_vm_double_cosine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sin(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_sine(A.size, &A[0], &res[0]) + clib.test_vm_double_sine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_cos_I_sin(np.ndarray[double] A): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_cosine_I_sine(A.size, &A[0], &res[0]) + clib.test_vm_cosine_I_sine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_exp(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_exp(A.size, &A[0], &res[0], 1) + clib.test_vm_double_exp(int(A.size), &A[0], &res[0], 1) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_I_exp(np.ndarray[complex] A): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_exp(A.size, &A[0], &res[0]) + clib.test_vm_double_complex_exp(int(A.size), &A[0], &res[0]) return res From 1ffcebc8e6a1165239c498abb5f16647e6d78c3f Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Sun, 25 Aug 2024 17:39:34 -0400 Subject: [PATCH 17/27] test --- src/c_lib/base/base_model.pyx | 3 ++- src/c_lib/include/histogram.h | 16 ++++++++-------- src/c_lib/lib/histogram.c | 18 ++++++++---------- src/c_lib/lib/simulation.c | 23 +++++++++++++++++++++++ 4 files changed, 41 insertions(+), 19 deletions(-) diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index f2edb9140..e21a92fe6 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -191,6 +191,7 @@ def core_simulator(method, # numpy uint8 corresponds to an unsigned char C type # More types found here: https://numpy.org/doc/stable/user/basics.types.html#array-types-and-conversions-between-types cdef ndarray[cnp.uint8_t] f_contrib = np.asarray(freq_contrib, dtype=np.uint8) + print('f_contrib size', sizeof(f_contrib[0])) # affine transformation cdef ndarray[double] affine_matrix_c @@ -206,7 +207,7 @@ def core_simulator(method, affine_matrix_c[3] -= affine_matrix_c[1]*affine_matrix_c[2] # sites _______________________________________________________________________________ - cdef int number_of_sites, number_of_couplings + cdef unsigned int number_of_sites, number_of_couplings cdef ndarray[int] spin_index_ij cdef ndarray[float] spin_i cdef ndarray[double] gyromagnetic_ratio_i diff --git a/src/c_lib/include/histogram.h b/src/c_lib/include/histogram.h index 6eb57029f..f759a216c 100644 --- a/src/c_lib/include/histogram.h +++ b/src/c_lib/include/histogram.h @@ -9,19 +9,19 @@ extern void histogram1d_c(int x_count, const double x_min, const double x_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict weights, const int stride_w); extern void histogram2d_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w); extern void histogram2d_interp_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w); diff --git a/src/c_lib/lib/histogram.c b/src/c_lib/lib/histogram.c index 380c8be54..9cefc1e9e 100644 --- a/src/c_lib/lib/histogram.c +++ b/src/c_lib/lib/histogram.c @@ -13,9 +13,8 @@ #include void histogram1d_c(int x_count, const double x_min, const double x_max, - double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict weights, const int stride_w) { + double *restrict hist, int sample_count, double *restrict sample_x, + const int stride_x, double *restrict weights, const int stride_w) { int ix, counter = sample_count; double temp_sample, normx; double *x_pt = sample_x, *w_pt = weights; @@ -35,10 +34,9 @@ void histogram1d_c(int x_count, const double x_min, const double x_max, void histogram2d_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, - int sample_count, const double *restrict sample_x, - const int stride_x, const double *restrict sample_y, - const int stride_y, const double *restrict weights, - const int stride_w) { + int sample_count, double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w) { int ix, iy, hist_index, counter = sample_count; double temp_sample_x, temp_sample_y, normx, normy; double *x_pt = sample_x, *y_pt = sample_y, *w_pt = weights; @@ -66,9 +64,9 @@ void histogram2d_c(int x_count, const double x_min, const double x_max, int y_co void histogram2d_interp_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w) { + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w) { int ix, iy, p, counter = sample_count; bool left; double temp_sample_x, temp_sample_y, normx, normy; diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index 6159bc83d..6a1a1e896 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -59,6 +59,29 @@ void __mrsimulator_core( et al. `Computation of Orientational Averages in Solid-State NMR by Gaussian Spherical Quadrature` JMR, 132, 1998. https://doi.org/10.1006/jmre.1998.1427 */ + printf("int %lu\n", sizeof(int)); + printf("double %lu\n", sizeof(double)); + printf("float %lu\n", sizeof(float)); + printf("unsigned int %lu\n", sizeof(unsigned int)); + printf("unsigned char %lu\n\n", sizeof(unsigned char)); + + printf("transition_pathway (float) %lu\n", sizeof(transition_pathway)); + printf("transition_pathway_weight (double) %lu\n", sizeof(transition_pathway_weight)); + printf("n_dimension (int) %lu\n\n", sizeof(n_dimension)); + + printf("iso_intrp (unsigned int) %lu\n", sizeof(iso_intrp)); + printf("freq_contrib (unsigned char) %lu\n", sizeof(freq_contrib)); + printf("affine_matrix double %lu\n\n", sizeof(affine_matrix)); + + printf("n_sites site (unsigned int) %lu\n", sizeof(sites->number_of_sites)); + printf("n_sites spin (float) %lu\n", sizeof(sites->spin)); + printf("n_sites gyromagnetic_ratio (double) %lu\n\n", + sizeof(sites->gyromagnetic_ratio)); + + printf("n_coupling n_c (unsigned int) %lu\n", sizeof(couplings->number_of_couplings)); + printf("n_coupling site_index (int) %lu\n", sizeof(couplings->site_index)); + printf("isotropic_j_in_Hz (double) %lu\n\n", sizeof(couplings->isotropic_j_in_Hz)); + unsigned char is_spectral; unsigned int evt; int dim, total_pts = scheme->n_gamma * scheme->total_orientations; From 76eed11d002f38ca0b81a43a697d4ab1570022b5 Mon Sep 17 00:00:00 2001 From: unknown <21365911+deepanshs@users.noreply.github.com> Date: Sun, 25 Aug 2024 20:42:38 -0400 Subject: [PATCH 18/27] test --- src/c_lib/base/base_model.pyx | 71 ++++++++++++------ src/c_lib/lib/simulation.c | 46 ++++++------ src/mrsimulator/simulator/__init__.py | 73 +++++++++++-------- .../test_lineshapes.py | 58 ++++++++++----- 4 files changed, 153 insertions(+), 95 deletions(-) diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index e21a92fe6..f90944c94 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -191,7 +191,6 @@ def core_simulator(method, # numpy uint8 corresponds to an unsigned char C type # More types found here: https://numpy.org/doc/stable/user/basics.types.html#array-types-and-conversions-between-types cdef ndarray[cnp.uint8_t] f_contrib = np.asarray(freq_contrib, dtype=np.uint8) - print('f_contrib size', sizeof(f_contrib[0])) # affine transformation cdef ndarray[double] affine_matrix_c @@ -207,7 +206,7 @@ def core_simulator(method, affine_matrix_c[3] -= affine_matrix_c[1]*affine_matrix_c[2] # sites _______________________________________________________________________________ - cdef unsigned int number_of_sites, number_of_couplings + cdef unsigned int number_of_sites, number_of_couplings, _site_idx_, _coup_idx_ cdef ndarray[int] spin_index_ij cdef ndarray[float] spin_i cdef ndarray[double] gyromagnetic_ratio_i @@ -279,22 +278,22 @@ def core_simulator(method, # Extract and assign site information from Site objects to C structure # --------------------------------------------------------------------- - for i in range(number_of_sites): - site = spin_sys.sites[i] - spin_i[i] = site.isotope.spin - gyromagnetic_ratio_i[i] = site.isotope.gyromagnetic_ratio - i3 = 3*i + for _site_idx_ in range(number_of_sites): + site = spin_sys.sites[_site_idx_] + spin_i[_site_idx_] = np.float32(site.isotope.spin) + gyromagnetic_ratio_i[_site_idx_] = site.isotope.gyromagnetic_ratio + i3 = 3 * _site_idx_ # CSA tensor if site.isotropic_chemical_shift is not None: - iso_n[i] = site.isotropic_chemical_shift + iso_n[_site_idx_] = site.isotropic_chemical_shift shielding = site.shielding_symmetric if shielding is not None: if shielding.zeta is not None: - zeta_n[i] = shielding.zeta + zeta_n[_site_idx_] = shielding.zeta if shielding.eta is not None: - eta_n[i] = shielding.eta + eta_n[_site_idx_] = shielding.eta if shielding.alpha is not None: ori_n[i3] = shielding.alpha if shielding.beta is not None: @@ -312,9 +311,9 @@ def core_simulator(method, quad = site.quadrupolar if quad is not None: if quad.Cq is not None: - Cq_e[i] = quad.Cq + Cq_e[_site_idx_] = quad.Cq if quad.eta is not None: - eta_e[i] = quad.eta + eta_e[_site_idx_] = quad.eta if quad.alpha is not None: ori_e[i3] = quad.alpha if quad.beta is not None: @@ -360,21 +359,21 @@ def core_simulator(method, ori_d = np.zeros(3*number_of_couplings, dtype=np.float64) # Extract and assign coupling information from Site objects to C structure - for i in range(number_of_couplings): - coupling = spin_sys.couplings[i] - spin_index_ij[2*i: 2*i+2] = coupling.site_index - i3 = 3*i + for _coup_idx_ in range(number_of_couplings): + coupling = spin_sys.couplings[_coup_idx_] + spin_index_ij[2*_coup_idx_: 2*_coup_idx_+2] = coupling.site_index + i3 = 3 * _coup_idx_ # J tensor if coupling.isotropic_j is not None: - iso_j[i] = coupling.isotropic_j + iso_j[_coup_idx_] = coupling.isotropic_j J_sym = coupling.j_symmetric if J_sym is not None: if J_sym.zeta is not None: - zeta_j[i] = J_sym.zeta + zeta_j[_coup_idx_] = J_sym.zeta if J_sym.eta is not None: - eta_j[i] = J_sym.eta + eta_j[_coup_idx_] = J_sym.eta if J_sym.alpha is not None: ori_j[i3] = J_sym.alpha if J_sym.beta is not None: @@ -386,9 +385,9 @@ def core_simulator(method, dipolar = coupling.dipolar if dipolar is not None: if dipolar.D is not None: - D_d[i] = dipolar.D + D_d[_coup_idx_] = dipolar.D if dipolar.eta is not None: - eta_d[i] = dipolar.eta + eta_d[_coup_idx_] = dipolar.eta if dipolar.alpha is not None: ori_d[i3] = dipolar.alpha if dipolar.beta is not None: @@ -448,6 +447,35 @@ def core_simulator(method, # print('pathway', transition_pathway_c) # print('weight', transition_pathway_weight_c) # print('pathway_count, inc', pathway_count, pathway_increment) + + if debug: + print("int ", sizeof(int)); + print("double ", sizeof(double)); + print("float ", sizeof(float)); + print("unsigned int ", sizeof(unsigned int)); + print("unsigned char ", sizeof(unsigned char)); + print() + + print("transition_pathway (float) ", sizeof(transition_pathway_c[0])); + print("transition_pathway_weight (double) ", sizeof(transition_pathway_weight_c[0])); + print("n_dimension (int) ", sizeof(n_dimension)); + print() + + print("iso_intrp (unsigned int) ", sizeof(isotropic_interpolation)); + print("freq_contrib (unsigned char) ", sizeof(f_contrib[0])); + print("affine_matrix (double) ", sizeof(affine_matrix_c[0])); + print() + + print("n_sites site (unsigned int) ", sizeof(sites_c.number_of_sites)); + print("n_sites spin (float) ", sizeof(sites_c.spin[0])); + print("n_sites gyromagnetic_ratio (double) ", sizeof(sites_c.gyromagnetic_ratio[0])); + print() + + print("n_coupling n_c (unsigned int) ", sizeof(couplings_c.number_of_couplings)); + print("n_coupling site_index (int) ", sizeof(couplings_c.site_index[0])); + print("isotropic_j_in_Hz (double) ", sizeof(couplings_c.isotropic_j_in_Hz[0])); + print() + for trans__ in range(pathway_count): clib.__mrsimulator_core( &[0], # as complex array @@ -483,6 +511,7 @@ def core_simulator(method, amp1.shape = method.shape() if gyromagnetic_ratio < 0: amp1 = np.fft.fftn(np.fft.ifftn(amp1).conj()) + amp1 = [amp1] clib.MRS_free_dimension(dimensions, n_dimension) clib.MRS_free_averaging_scheme(averaging_scheme) diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index 6a1a1e896..11f96e6aa 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -59,28 +59,30 @@ void __mrsimulator_core( et al. `Computation of Orientational Averages in Solid-State NMR by Gaussian Spherical Quadrature` JMR, 132, 1998. https://doi.org/10.1006/jmre.1998.1427 */ - printf("int %lu\n", sizeof(int)); - printf("double %lu\n", sizeof(double)); - printf("float %lu\n", sizeof(float)); - printf("unsigned int %lu\n", sizeof(unsigned int)); - printf("unsigned char %lu\n\n", sizeof(unsigned char)); - - printf("transition_pathway (float) %lu\n", sizeof(transition_pathway)); - printf("transition_pathway_weight (double) %lu\n", sizeof(transition_pathway_weight)); - printf("n_dimension (int) %lu\n\n", sizeof(n_dimension)); - - printf("iso_intrp (unsigned int) %lu\n", sizeof(iso_intrp)); - printf("freq_contrib (unsigned char) %lu\n", sizeof(freq_contrib)); - printf("affine_matrix double %lu\n\n", sizeof(affine_matrix)); - - printf("n_sites site (unsigned int) %lu\n", sizeof(sites->number_of_sites)); - printf("n_sites spin (float) %lu\n", sizeof(sites->spin)); - printf("n_sites gyromagnetic_ratio (double) %lu\n\n", - sizeof(sites->gyromagnetic_ratio)); - - printf("n_coupling n_c (unsigned int) %lu\n", sizeof(couplings->number_of_couplings)); - printf("n_coupling site_index (int) %lu\n", sizeof(couplings->site_index)); - printf("isotropic_j_in_Hz (double) %lu\n\n", sizeof(couplings->isotropic_j_in_Hz)); + // printf("int %zu\n", sizeof(int)); + // printf("double %zu\n", sizeof(double)); + // printf("float %zu\n", sizeof(float)); + // printf("unsigned int %zu\n", sizeof(unsigned int)); + // printf("unsigned char %zu\n\n", sizeof(unsigned char)); + + // printf("transition_pathway (float) %zu\n", sizeof(transition_pathway[0])); + // printf("transition_pathway_weight (double) %zu\n", + // sizeof(transition_pathway_weight[0])); printf("n_dimension (int) %zu\n\n", + // sizeof(n_dimension)); + + // printf("iso_intrp (unsigned int) %zu\n", sizeof(iso_intrp)); + // printf("freq_contrib (unsigned char) %zu\n", sizeof(freq_contrib[0])); + // printf("affine_matrix double %zu\n\n", sizeof(affine_matrix[0])); + + // printf("n_sites site (unsigned int) %zu\n", sizeof(sites->number_of_sites)); + // printf("n_sites spin (float) %zu\n", sizeof(sites->spin[0])); + // printf("n_sites gyromagnetic_ratio (double) %zu\n\n", + // sizeof(sites->gyromagnetic_ratio[0])); + + // printf("n_coupling n_c (unsigned int) %zu\n", + // sizeof(couplings->number_of_couplings)); printf("n_coupling site_index (int) + // %zu\n", sizeof(couplings->site_index[0])); printf("isotropic_j_in_Hz (double) + // %zu\n\n", sizeof(couplings->isotropic_j_in_Hz[0])); unsigned char is_spectral; unsigned int evt; diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 3f1e9a6d4..82ff585d2 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -8,8 +8,6 @@ import numpy as np import pandas as pd import psutil -from joblib import delayed -from joblib import Parallel from mrsimulator import Site from mrsimulator import SpinSystem from mrsimulator.base_model import core_simulator @@ -22,6 +20,9 @@ from .config import ConfigSimulator +# from joblib import delayed +# from joblib import Parallel + __author__ = "Deepansh Srivastava" __email__ = "srivastava.89@osu.edu" @@ -383,7 +384,7 @@ def run( >>> sim.run() # doctest:+SKIP """ - verbose = 0 + # verbose = 0 if opt is None: opt = self.optimize() @@ -395,29 +396,39 @@ def run( for index in method_index: method = self.methods[index] - spin_sys = get_chunks(self.spin_systems, n_jobs) - - # Chunking transition pathways and weights form the optimization dictionary - pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) - weights = get_chunks(opt["precomputed_weights"][index], n_jobs) config_dict = self.config.get_int_dict() - jobs = ( - delayed(core_simulator)( - method=method, - spin_systems=sys, - transition_pathways=pth, - transition_weights=wht, - **config_dict, - **kwargs, - ) - for sys, pth, wht in zip(spin_sys, pathways, weights) + amp = core_simulator( + method=method, + spin_systems=self.spin_systems, + transition_pathways=opt["precomputed_pathways"][index], + transition_weights=opt["precomputed_weights"][index], + **config_dict, + **kwargs, ) - amp = Parallel( - n_jobs=n_jobs, - verbose=verbose, - backend="loky", - )(jobs) + + # spin_sys = get_chunks(self.spin_systems, n_jobs) + + # # Chunk transition pathways and weights form the optimization dictionary + # pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) + # weights = get_chunks(opt["precomputed_weights"][index], n_jobs) + + # jobs = ( + # delayed(core_simulator)( + # method=method, + # spin_systems=sys, + # transition_pathways=pth, + # transition_weights=wht, + # **config_dict, + # **kwargs, + # ) + # for sys, pth, wht in zip(spin_sys, pathways, weights) + # ) + # amp = Parallel( + # n_jobs=n_jobs, + # verbose=verbose, + # backend="loky", + # )(jobs) gyromagnetic_ratio = method.channels[0].gyromagnetic_ratio B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density @@ -575,17 +586,15 @@ def _package_amp_after_simulation(self, method, amp, pack_as_csdm): bool pack_as_csdm: Packages the simulated spectrum as a CSDM object if true. Otherwise kept as a numpy array. """ - if isinstance(amp[0], list): - simulated_dataset = [] - for item in amp: - simulated_dataset += item - if isinstance(amp[0], np.ndarray): - simulated_dataset = [np.asarray(amp).sum(axis=0)] + # if isinstance(amp[0], list): + # simulated_dataset = [] + # for item in amp: + # simulated_dataset += item + # if isinstance(amp[0], np.ndarray): + # simulated_dataset = [np.asarray(amp).sum(axis=0)] method.simulation = ( - self._as_csdm_object(simulated_dataset, method) - if pack_as_csdm - else np.asarray(simulated_dataset) + self._as_csdm_object(amp, method) if pack_as_csdm else np.asarray(amp) ) diff --git a/tests/spectral_integration_tests/test_lineshapes.py b/tests/spectral_integration_tests/test_lineshapes.py index f22a2117f..46e6158c2 100644 --- a/tests/spectral_integration_tests/test_lineshapes.py +++ b/tests/spectral_integration_tests/test_lineshapes.py @@ -149,7 +149,8 @@ def test_pure_shielding_sideband_simpson(report): "failed to compare shielding sidebands with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "shielding_sidebands") - for i in range(8): + max_l = 8 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] @@ -163,7 +164,9 @@ def test_pure_shielding_sideband_simpson(report): ) res.append([data_mrsimulator, data_source]) - compile_plots(dim, res, info, title="Shielding Sidebands", report=report) + compile_plots( + dim, res, info, title=f"Shielding Sidebands {i}/{max_l}", report=report + ) message = f"{error_message} test0{i}.json" check_all_close(res, message, rel_limit=1e-3) @@ -174,7 +177,8 @@ def test_pure_quadrupolar_sidebands_simpson(report): "failed to compare quadrupolar sidebands with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "quad_sidebands") - for i in range(2): + max_l = 2 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") for volume in VOLUMES: res = [] @@ -189,7 +193,9 @@ def test_pure_quadrupolar_sidebands_simpson(report): ) res.append([data_mrsimulator, data_source]) - compile_plots(dim, res, info, title="Quad Sidebands", report=report) + compile_plots( + dim, res, info, title=f"Quad Sidebands {i}/{max_l}", report=report + ) message = f"{error_message} test0{i:02d}.json" check_all_close(res, message, rel_limit=1e-3) @@ -200,7 +206,8 @@ def test_csa_plus_quadrupolar_lineshape_simpson(report): "failed to compare quad + csa lineshape with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "csa_quad") - for i in range(10): + max_l = 10 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] for volume in VOLUMES: @@ -213,7 +220,7 @@ def test_csa_plus_quadrupolar_lineshape_simpson(report): dim, res, info, - title="Quad + Shielding Sidebands", + title=f"Quad + Shielding Sidebands {i}/{max_l}", report=report, ) @@ -226,7 +233,8 @@ def test_1st_order_quadrupolar_lineshape_simpson(report): "failed to compare 1st order quad lineshape with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "quad_1st_order") - for i in range(2): + max_l = 2 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] for volume in VOLUMES: @@ -239,7 +247,7 @@ def test_1st_order_quadrupolar_lineshape_simpson(report): dim, res, info, - title="1st Order Quadrupolar Lineshape", + title=f"1st Order Quadrupolar Lineshape {i}/{max_l}", report=report, ) @@ -250,7 +258,8 @@ def test_1st_order_quadrupolar_lineshape_simpson(report): def test_j_coupling_lineshape_simpson(report): error_message = "failed to compare j-coupling with simpson simulation from file" path_ = path.join(SIMPSON_TEST_PATH, "j-coupling") - for i in range(20): + max_l = 20 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] for volume in VOLUMES: @@ -259,7 +268,9 @@ def test_j_coupling_lineshape_simpson(report): ) res.append([data_mrsimulator, data_source]) - compile_plots(dim, res, info, title="J-coupling Spectra", report=report) + compile_plots( + dim, res, info, title=f"J-coupling Spectra {i}/{max_l}", report=report + ) message = f"{error_message} test0{i:02d}.json" check_all_close(res, message, rel_limit=9e-3) @@ -270,7 +281,8 @@ def test_dipolar_coupling_lineshape_simpson(report): "failed to compare dipolar-coupling with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "dipolar-coupling") - for i in range(7): + max_l = 7 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] for volume in VOLUMES: @@ -279,7 +291,9 @@ def test_dipolar_coupling_lineshape_simpson(report): ) res.append([data_mrsimulator, data_source]) - compile_plots(dim, res, info, title="Dipolar-coupling Spectra", report=report) + compile_plots( + dim, res, info, title=f"Dipolar-coupling Spectra {i}/{max_l}", report=report + ) message = f"{error_message} test0{i:02d}.json" check_all_close(res, message, rel_limit=4e-3) @@ -290,7 +304,8 @@ def test_2D_sideband_sideband_simpson(report): "failed to compare sideband-sideband with simpson simulation from file" ) path_ = path.join(SIMPSON_TEST_PATH, "sideband_sideband") - for i in range(5): + max_l = 5 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] @@ -303,7 +318,7 @@ def test_2D_sideband_sideband_simpson(report): res, info, dim2=dim[1], - title="2D sideband-sideband", + title=f"2D sideband-sideband {i}/{max_l}", report=report, label="Simpson", ) @@ -321,9 +336,10 @@ def test_quad_csa_cross_rmnsim(report): error_message = ( "failed to compare quad-csa cross-term spectra with rmnsim from file" ) + max_l = 7 for folder in ["quad_csa_cross1", "quad_csa_cross2"]: path_ = path.join(RNMSIM_TEST_PATH, folder) - for i in range(7): + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] @@ -337,7 +353,7 @@ def test_quad_csa_cross_rmnsim(report): dim, res, info, - title="Quad-CSA 2nd Order Cross-Term", + title=f"Quad-CSA 2nd Order Cross-Term{i}/{max_l}", report=report, label="rmnsim", ) @@ -357,7 +373,8 @@ def test_pure_shielding_static_lineshape_python_brute(report): "failed to compare shielding lineshape with brute force simulation from file" ) path_ = path.join(PYTHON_BRUTE_TEST_PATH, "shielding") - for i in range(5): + max_l = 5 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] @@ -375,7 +392,7 @@ def test_pure_shielding_static_lineshape_python_brute(report): dim, res, info, - title="Shielding Static Lineshape (Brute Force)", + title=f"Shielding Static Lineshape (Brute Force) {i}/{max_l}", report=report, label="Brute", ) @@ -394,7 +411,8 @@ def test_pure_quadrupolar_lineshape_python_brute(report): "failed to compare quadrupolar lineshape with brute force simulation from file" ) path_ = path.join(PYTHON_BRUTE_TEST_PATH, "quad") - for i in range(19): + max_l = 19 + for i in range(max_l): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") res = [] @@ -411,7 +429,7 @@ def test_pure_quadrupolar_lineshape_python_brute(report): dim, res, info, - title="Quad Lineshape Self-Test", + title=f"Quad Lineshape Self-Test {i}/{max_l}", report=report, label="self", ) From 24d087ff12c796a4d9c5aa893aebec603e42c046 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Tue, 27 Aug 2024 06:07:10 -0400 Subject: [PATCH 19/27] clean interpolation code --- src/c_lib/include/interpolation.h | 47 ++-- src/c_lib/include/vm.h | 5 +- src/c_lib/lib/frequency_averaging.c | 4 +- src/c_lib/lib/interpolation.c | 207 +++++++++--------- src/c_lib/test/test.pxd | 44 ++-- src/c_lib/test/test.pyx | 34 ++- ...orientations_and_triangle_interpolation.py | 3 +- 7 files changed, 168 insertions(+), 176 deletions(-) diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index 1f8d2cc94..bde407a74 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -12,21 +12,21 @@ /** * @brief Create a triangle with coordinates (f1, f2, f2) onto a 1D grid. * - * @param f1 A pointer to the coordinate f11. - * @param f2 A pointer to the coordinate f12. - * @param f3 A pointer to the coordinate f13. - * @param amp A pointer to the area of the vector. + * @param f1 Coordinate f11. + * @param f2 Coordinate f12. + * @param f3 Coordinate f13. + * @param amp Area of triangle. * @param spec A pointer to the starting index of a one-dimensional array. - * @param m0 A pointer to the number of points on the 1D grid. + * @param m0 The number of points on the 1D grid. * @param type The type of interpolation for the isotropic components. * 0. delta-interpolation, * 1. gaussian-interpolation. */ -extern void triangle_interpolation1D(double *f1, double *f2, double *f3, double *amp, - double *spec, int *m0, unsigned int iso_intrp); +extern void triangle_interpolation1D(double f1, double f2, double f3, double amp, + double *spec, int m0, unsigned int iso_intrp); -extern void triangle_interpolation1D_linear(double *f1, double *f2, double *f3, - double *amp, double *spec, int *m0); +extern void triangle_interpolation1D_linear(double f1, double f2, double f3, double amp, + double *spec, int m0); extern void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, double *amp_real, double *amp_imag, int dimension_count, @@ -60,42 +60,41 @@ extern void generic_1d_triangle_average(double *spec, const unsigned int freq_si const unsigned int position_size, int32_t *positions, const unsigned int nt); -extern void triangle_interpolation1D_gaussian(double *f1, double *f2, double *f3, - double *amp, double *spec, int *m0); +extern void triangle_interpolation1D_gaussian(double f1, double f2, double f3, + double amp, double *spec, int m0); /** * @brief Rasterize a vector triangle with coordinates ((f11, f21), (f12, f22), (f13, * f23)) onto a 2D grid. * - * @param f11 A pointer to the coordinate f11. - * @param f12 A pointer to the coordinate f12. - * @param f13 A pointer to the coordinate f13. - * @param f21 A pointer to the coordinate f21. - * @param f22 A pointer to the coordinate f22. - * @param f23 A pointer to the coordinate f23. - * @param amp A pointer to the area of the vector. + * @param f11 Coordinate f11. + * @param f12 Coordinate f12. + * @param f13 Coordinate f13. + * @param f21 Coordinate f21. + * @param f22 Coordinate f22. + * @param f23 Coordinate f23. + * @param amp Area of 2D triangle. * @param spec A pointer to the starting index of a two-dimensional array. * @param m0 An integer with the rows in the 2D grid. * @param m1 An integer with the columns in the 2D grid. * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ -extern void triangle_interpolation2D(double *f11, double *f12, double *f13, double *f21, - double *f22, double *f23, double *amp, - double *spec, int m0, int m1, - unsigned int iso_intrp); +extern void triangle_interpolation2D(double f11, double f12, double f13, double f21, + double f22, double f23, double amp, double *spec, + int m0, int m1, unsigned int iso_intrp); /** * @brief Sum amplitudes from the triangles interpolations over the region of an octant. * The samplings over the octant is as per Alderman and Grand scheme. * * @param nt Number of triangles along the edge of the octant. - * @param freq A pointer to an array of frequencies evaluated at octant coordinates. + * @param freq Delta frequency. * @param amp A pointer to the amplitudes for the frequencies at octant coordinates. * @param stride Stride setp for the amplitudes (amp) array. * @param n_spec Number of points in the spectrum array (spec) * @param spec A pointer to the starting index of a one-dimensional array * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ -void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *amp, +void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, int stride, int n_spec, double *spec, unsigned int iso_intrp); diff --git a/src/c_lib/include/vm.h b/src/c_lib/include/vm.h index 020c477d1..b4156a17a 100644 --- a/src/c_lib/include/vm.h +++ b/src/c_lib/include/vm.h @@ -217,8 +217,7 @@ static inline void vm_double_square_root_inplace(int count, double *restrict x) * res = x * y */ static inline void vm_double_complex_multiply(int count, const void *restrict x, - const void *restrict y, - void *restrict res) { + const void *restrict y, void *res) { double *res_ = (double *)res; double *x_ = (double *)x; double *y_ = (double *)y; @@ -440,7 +439,7 @@ static inline void vm_double_complex_exp(int count, const void *restrict x, * res = exp(x(imag)) */ static inline void vm_double_complex_exp_imag_only(int count, const void *restrict x, - void *restrict res) { + void *res) { // x = __builtin_assume_aligned(x, 32); // res = __builtin_assume_aligned(res, 32); double *x_ = (double *)x; diff --git a/src/c_lib/lib/frequency_averaging.c b/src/c_lib/lib/frequency_averaging.c index fc2c5557f..09a307e38 100644 --- a/src/c_lib/lib/frequency_averaging.c +++ b/src/c_lib/lib/frequency_averaging.c @@ -120,9 +120,9 @@ void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * vm_double_multiply(scheme->total_orientations, phase_ptr + 1, 2, &s[k1], amps_imag); while (j++ < planA->n_octants) { - octahedronDeltaInterpolation(nt, &offset, &s_real[address], 1, + octahedronDeltaInterpolation(nt, offset, &s_real[address], 1, dimensions->count, spec, iso_intrp); - octahedronDeltaInterpolation(nt, &offset, &s_imag[address], 1, + octahedronDeltaInterpolation(nt, offset, &s_imag[address], 1, dimensions->count, spec + 1, iso_intrp); address += npts; } diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index 15c31a984..ac807bbe8 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -23,51 +23,51 @@ * @param amp The area corresponding to the frequency coordinate. * @param spec1 A pointer to the intensity vector. */ -static void inline delta_fn_linear_interpolation(const double *freq, const int *points, - double *amp, double *spec1) { - int p = (int)*freq; - if (p >= *points || p < 0) return; +static void inline delta_fn_linear_interpolation(const double freq, const int points, + double amp, double *spec1) { + int p = (int)freq; + if (p >= points || p < 0) return; double diff, delta; bool left; - diff = *freq - (double)p; + diff = freq - (double)p; delta = diff - 0.5; double *spec = &spec1[2 * p]; // Do not interpolate the intensity if the difference < 1.0e-6. // Ensures that the sideband dimension frequencies do not get interpolated. if (absd(delta) < TOL) { - *spec += *amp; + *spec += amp; return; } // Linear interpolation left = delta < 0; if (left) { - if (p != 0) *(spec - 2) -= *amp * delta; - *spec += *amp * (1.0 + delta); + if (p != 0) *(spec - 2) -= amp * delta; + *spec += amp * (1.0 + delta); } else { - if (p + 1 != *points) *(spec + 2) += *amp * delta; - *spec += *amp * (1.0 - delta); + if (p + 1 != points) *(spec + 2) += amp * delta; + *spec += amp * (1.0 - delta); } } -static void inline delta_fn_gauss_interpolation(const double *freq, const int *points, - double *amp, double *spec) { +static void inline delta_fn_gauss_interpolation(const double freq, const int points, + double amp, double *spec) { double res, w, a0, a1, a2, a3, a4, sum, temp; - int p = (int)(floor(*freq)), index, ip2, ip1, im1, im2, pad = 2, p2 = 2 * p; - if (p >= *points + pad || p < -pad + 1) return; + int p = (int)(floor(freq)), index, ip2, ip1, im1, im2, pad = 2, p2 = 2 * p; + if (p >= points + pad || p < -pad + 1) return; // for sideband delta freq. It avoids round-off errors. - res = *freq - (double)p; - if (absd(res - 0.5) < TOL && p >= 0 && p < *points) { - spec[p2] += *amp; + res = freq - (double)p; + if (absd(res - 0.5) < TOL && p >= 0 && p < points) { + spec[p2] += amp; return; } - p = (int)(floor(*freq - 0.5)); + p = (int)(floor(freq - 0.5)); p2 = 2 * p; - res = *freq - (double)p - 0.5; + res = freq - (double)p - 0.5; res *= gauss_table_precision_inverse; index = (int)res; w = res - (double)index; @@ -89,18 +89,18 @@ static void inline delta_fn_gauss_interpolation(const double *freq, const int *p sum += a3; sum += a4; - temp = *amp / sum; + temp = amp / sum; // double *spec = &spec1[2 * p]; - if (p - 2 >= 0 && p - 2 < *points) spec[p2 - 4] += temp * a0; - if (p - 1 >= 0 && p - 1 < *points) spec[p2 - 2] += temp * a1; - if (p >= 0 && p < *points) spec[p2] += temp * a2; - if (p + 1 >= 0 && p + 1 < *points) spec[p2 + 2] += temp * a3; - if (p + 2 >= 0 && p + 2 < *points) spec[p2 + 4] += temp * a4; + if (p - 2 >= 0 && p - 2 < points) spec[p2 - 4] += temp * a0; + if (p - 1 >= 0 && p - 1 < points) spec[p2 - 2] += temp * a1; + if (p >= 0 && p < points) spec[p2] += temp * a2; + if (p + 1 >= 0 && p + 1 < points) spec[p2 + 2] += temp * a3; + if (p + 2 >= 0 && p + 2 < points) spec[p2 + 4] += temp * a4; } /** - * @brief Get the clipping conditions object + * @brief Get the clipping conditions object. Update the value of p, pmid, and pmax * * @param p Start index of the triangle. * @param pmid Apex index of the triangle. @@ -108,19 +108,19 @@ static void inline delta_fn_gauss_interpolation(const double *freq, const int *p * @param points Max attainable index. * @param clips Pointer to bool array to store the clip conditions. */ -static inline void get_clipping_conditions(int *p, int *pmid, int *pmax, int *points, +static inline void get_clipping_conditions(int *p, int *pmid, int *pmax, int points, bool *clips) { *clips = *p <= 0; // left triangle left clip *p = *clips++ ? 0 : *p; - *clips = *pmid > *points - 1; // left triangle right clip - *pmid = *clips++ ? *points - 1 : *pmid; + *clips = *pmid > points - 1; // left triangle right clip + *pmid = *clips++ ? points - 1 : *pmid; *clips = *pmid <= 0; // right triangle left clip *pmid = *clips++ ? 0 : *pmid; - *clips = *pmax > *points - 1; // right triangle right clip - *pmax = *clips ? *points - 1 : *pmax; + *clips = *pmax > points - 1; // right triangle right clip + *pmax = *clips ? points - 1 : *pmax; } /** @@ -210,20 +210,20 @@ static inline void right_triangle_interpolate(int p, int pmax, bool l_clip, bool *(spec += 2) += (!r_clip) ? df * df * 0.5 * df1 : diff - df1; } -static inline void __triangle_interpolation(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { +static inline void __triangle_interpolation(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { int p, pmid, pmax, i, j; double top, t; - p = (int)(*freq1); + p = (int)(freq1); // check if the three points lie within a bin interval. - if (p == (int)(*freq2) && p == (int)(*freq3)) { - if (p >= *points || p < 0) return; - spec[2 * p] += *amp; + if (p == (int)(freq2) && p == (int)(freq3)) { + if (p >= points || p < 0) return; + spec[2 * p] += amp; return; } - double f[3] = {freq1[0], freq2[0], freq3[0]}; + double f[3] = {freq1, freq2, freq3}; // arrange the numbers in ascending order (sort) @@ -239,14 +239,14 @@ static inline void __triangle_interpolation(double *freq1, double *freq2, double // if min frequency is higher than the last bin, return p = (int)f[0]; - if (p >= *points) return; + if (p >= points) return; // if max frequency is lower than the first bin, return pmax = (int)f[2]; if (pmax < 0) return; pmid = (int)f[1]; - top = *amp * 2.0 / (f[2] - f[0]); + top = amp * 2.0 / (f[2] - f[0]); bool clips[4]; get_clipping_conditions(&p, &pmid, &pmax, points, clips); @@ -260,9 +260,9 @@ static inline void __triangle_interpolation(double *freq1, double *freq2, double } } -void triangle_interpolation1D(double *freq1, double *freq2, double *freq3, double *amp, - double *spec, int *points, unsigned int iso_intrp) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { +void triangle_interpolation1D(double freq1, double freq2, double freq3, double amp, + double *spec, int points, unsigned int iso_intrp) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { if (iso_intrp == 0) delta_fn_linear_interpolation(freq1, points, amp, spec); if (iso_intrp == 1) delta_fn_gauss_interpolation(freq1, points, amp, spec); } else { @@ -270,18 +270,18 @@ void triangle_interpolation1D(double *freq1, double *freq2, double *freq3, doubl } } -void triangle_interpolation1D_linear(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { +void triangle_interpolation1D_linear(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { delta_fn_linear_interpolation(freq1, points, amp, spec); } else { __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); } } -void triangle_interpolation1D_gaussian(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { +void triangle_interpolation1D_gaussian(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { delta_fn_gauss_interpolation(freq1, points, amp, spec); } else { __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); @@ -355,7 +355,7 @@ void triangle_interpolation1D_gaussian(double *freq1, double *freq2, double *fre // / \ _____ \. // / \_____ \. // f10 ------------------ f11 bottom line -static inline void quadrilateral_bin(double *f00, double *f11, double *f10, double *f01, +static inline void quadrilateral_bin(double f00, double f11, double f10, double f01, double top, double bottom, double total, double amp, double *spec, int m1, unsigned int iso_intrp) { @@ -365,10 +365,10 @@ static inline void quadrilateral_bin(double *f00, double *f11, double *f10, doub norm = amp / total; amp_bottom = bottom * norm; amp_top = top * norm; - triangle_interpolation1D(f00, f11, f10, &_bottom, spec, &m1, iso_intrp); - triangle_interpolation1D(f00, f11, f01, &_top, spec, &m1, iso_intrp); + triangle_interpolation1D(f00, f11, f10, amp_bottom, spec, m1, iso_intrp); + triangle_interpolation1D(f00, f11, f01, amp_top, spec, m1, iso_intrp); } else { - triangle_interpolation1D(f00, f11, f10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f00, f11, f10, amp, spec, m1, iso_intrp); } } @@ -397,11 +397,11 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * f10 + f2[0]; if (!r_clip && !l_clip) { amp = f10 * top * 0.5; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } if (!r_clip_r && l_clip) { amp = f1[1] * (f10 - f1[0]) * 0.5 * df1; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } return; } @@ -416,7 +416,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * diff + f2[0]; if (!l_clip) { amp = 0.5 * diff * diff * df1; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } else { amp = (diff - 0.5) * df1; x00 = *x10 - f01_slope; @@ -424,7 +424,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, line_down = absd(*x11 - *x10); line_up = absd(x01 - x00); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } p++; @@ -444,7 +444,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, x01 = *x11; *x10 += f01_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); line_up += abs_sdiff; line_down += abs_sdiff; @@ -462,7 +462,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * f10 + f2[0]; line_down = absd(*x11 - *x10); denom = line_up + line_down; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } else { diff += df1; @@ -470,7 +470,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, x01 = *x11; *x10 += f01_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); } } @@ -499,11 +499,11 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, if (p == pmax) { if (!r_clip) { amp = f21 * top * 0.5; - triangle_interpolation1D(x11, &f2[1], &f2[2], &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x11, f2[1], f2[2], amp, spec, m1, iso_intrp); } if (r_clip && !r_clip_l) { amp = (f21 - diff) * (diff + f21) * 0.5 * df2; - triangle_interpolation1D(x11, &f2[1], &f2[2], &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x11, f2[1], f2[2], amp, spec, m1, iso_intrp); } return; } @@ -521,7 +521,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, line_up = absd(x01 - x00); line_down = absd(*x11 - *x10); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } else { amp = (diff + 0.5) * df2; @@ -532,7 +532,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, line_up = absd(x01 - x00); line_down = absd(*x11 - *x10); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } p++; @@ -550,7 +550,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, x01 = *x11; *x10 += f12_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); line_up -= abs_sdiff; line_down -= abs_sdiff; @@ -565,65 +565,64 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, amp = temp * temp * 0.5 * df2; x01 = *x11; *x11 = f2[2]; - triangle_interpolation1D(x10, x11, &x01, &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x10, *x11, x01, amp, spec, m1, iso_intrp); } else { diff -= df2; x00 = *x10; x01 = *x11; *x10 += f12_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); } } -void triangle_interpolation2D(double *freq11, double *freq12, double *freq13, - double *freq21, double *freq22, double *freq23, - double *amp, double *spec, int m0, int m1, - unsigned int iso_intrp) { +void triangle_interpolation2D(double freq11, double freq12, double freq13, + double freq21, double freq22, double freq23, double amp, + double *spec, int m0, int m1, unsigned int iso_intrp) { double top, t1, t2, diff, temp, n_i; int p, pmid, pmax, i, j; - double freq10_01, freq11_02; + double freq10_01 = 0.0, freq11_02 = 0.0; - p = (int)(*freq11); + p = (int)(freq11); - if (absd(freq11[0] - freq12[0]) < TOL && absd(freq11[0] - freq13[0]) < TOL) { + if (absd(freq11 - freq12) < TOL && absd(freq11 - freq13) < TOL) { if (p >= m0 || p < 0) return; - diff = freq11[0] - (double)p; + diff = freq11 - (double)p; n_i = 0.5; if (absd(diff - n_i) < TOL) { - triangle_interpolation1D(freq21, freq22, freq23, amp, &spec[2 * p * m1], &m1, + triangle_interpolation1D(freq21, freq22, freq23, amp, &spec[2 * p * m1], m1, iso_intrp); return; } if (diff < n_i) { if (p != 0) { - temp = *amp * (n_i - diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * (p - 1) * m1], - &m1, iso_intrp); + temp = amp * (n_i - diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * (p - 1) * m1], + m1, iso_intrp); } - temp = *amp * (n_i + diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * p * m1], &m1, + temp = amp * (n_i + diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * p * m1], m1, iso_intrp); return; } if (diff > n_i) { if (p + 1 != m0) { - temp = *amp * (diff - n_i); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * (p + 1) * m1], - &m1, iso_intrp); + temp = amp * (diff - n_i); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * (p + 1) * m1], + m1, iso_intrp); } - temp = *amp * (1 + n_i - diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * p * m1], &m1, + temp = amp * (1 + n_i - diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * p * m1], m1, iso_intrp); return; } return; } - double f1[3] = {freq11[0], freq12[0], freq13[0]}; - double f2[3] = {freq21[0], freq22[0], freq23[0]}; + double f1[3] = {freq11, freq12, freq13}; + double f2[3] = {freq21, freq22, freq23}; // arrange the numbers in ascending order for (j = 1; j <= 2; j++) { @@ -648,10 +647,10 @@ void triangle_interpolation2D(double *freq11, double *freq12, double *freq13, if (pmax < 0) return; pmid = (int)f1[1]; - top = *amp * 2.0 / (f1[2] - f1[0]); + top = amp * 2.0 / (f1[2] - f1[0]); bool clips[4]; - get_clipping_conditions(&p, &pmid, &pmax, &m0, clips); + get_clipping_conditions(&p, &pmid, &pmax, m0, clips); if (f1[1] >= 0.0) { lower_triangle_interpolation_2d(p, pmid, clips[0], clips[1], clips[2], top, f1, f2, &freq10_01, &freq11_02, m1, spec, iso_intrp); @@ -814,7 +813,7 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // self.x0 + t1*xdelta, self.y0 + t1*ydelta) // return [clipped_line] -void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *amp, +void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, int stride, int n_spec, double *spec, unsigned int iso_intrp) { int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; @@ -841,8 +840,8 @@ void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *a int_i_stride += stride; int_j_stride += stride; } - if (iso_intrp == 0) delta_fn_linear_interpolation(freq, &n_spec, &1, spec); - if (iso_intrp == 1) delta_fn_gauss_interpolation(freq, &n_spec, &1, spec); + if (iso_intrp == 0) delta_fn_linear_interpolation(freq, n_spec, amp1, spec); + if (iso_intrp == 1) delta_fn_gauss_interpolation(freq, n_spec, amp1, spec); } void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, @@ -862,12 +861,12 @@ void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, temp = amp[int_i_stride + stride] + amp_address[int_j_stride]; amp1 = temp + amp[int_i_stride]; - __triangle_interpolation(&freq[i], &freq[i + 1], &freq_address[j], &1, spec, &m); + __triangle_interpolation(freq[i], freq[i + 1], freq_address[j], amp1, spec, m); if (i < local_index) { temp += amp_address[int_j_stride + stride]; - __triangle_interpolation(&freq[i + 1], &freq_address[j], &freq_address[j + 1], - &temp, spec, &m); + __triangle_interpolation(freq[i + 1], freq_address[j], freq_address[j + 1], temp, + spec, m); } else { local_index = j + nt; i++; @@ -894,7 +893,7 @@ void generic_1d_triangle_interpolation(double *spec, const unsigned int freq_siz p3 = *positions++; // we do amp_sum because amps are already scaled to account for the factor 3 amp_sum = amp[p1] + amp[p2] + amp[p3]; - __triangle_interpolation(&freq[p1], &freq[p2], &freq[p3], &_sum, spec, &m); + __triangle_interpolation(freq[p1], freq[p2], freq[p3], amp_sum, spec, m); } } @@ -980,8 +979,8 @@ void generic_2d_triangle_interpolation(double *spec, const unsigned int freq_siz p3 = *positions++; // we do amp_sum because amps are already scaled to account for the factor 3 amp_sum = amp[amp_stride * p1] + amp[amp_stride * p2] + amp[amp_stride * p3]; - triangle_interpolation2D(&freq1[p1], &freq1[p2], &freq1[p3], &freq2[p1], &freq2[p2], - &freq2[p3], &_sum, spec, m0, m1, iso_intrp); + triangle_interpolation2D(freq1[p1], freq1[p2], freq1[p3], freq2[p1], freq2[p2], + freq2[p3], amp_sum, spec, m0, m1, iso_intrp); } } @@ -1035,15 +1034,15 @@ void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, int n temp = amp[int_i_stride + stride] + amp_address[int_j_stride]; amp1 = temp + amp[int_i_stride]; - triangle_interpolation2D(&freq1[i], &freq1[i + 1], &freq1_address[j], &freq2[i], - &freq2[i + 1], &freq2_address[j], &1, spec, m0, m1, + triangle_interpolation2D(freq1[i], freq1[i + 1], freq1_address[j], freq2[i], + freq2[i + 1], freq2_address[j], amp1, spec, m0, m1, iso_intrp); if (i < local_index) { temp += amp_address[int_j_stride + stride]; - triangle_interpolation2D(&freq1[i + 1], &freq1_address[j], &freq1_address[j + 1], - &freq2[i + 1], &freq2_address[j], &freq2_address[j + 1], - &temp, spec, m0, m1, iso_intrp); + triangle_interpolation2D(freq1[i + 1], freq1_address[j], freq1_address[j + 1], + freq2[i + 1], freq2_address[j], freq2_address[j + 1], + temp, spec, m0, m1, iso_intrp); } else { local_index = j + nt; i++; diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index a38629eb8..c7bda5af7 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -84,38 +84,38 @@ cdef extern from "octahedron.h": cdef extern from "interpolation.h": void triangle_interpolation1D( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points, + int points, unsigned int iso_intrp) void triangle_interpolation1D_linear( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points) + int points) void triangle_interpolation1D_gaussian( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points) + int points) void triangle_interpolation2D( - double *freq11, - double *freq12, - double *freq13, - double *freq21, - double *freq22, - double *freq23, - double *amp, + double freq11, + double freq12, + double freq13, + double freq21, + double freq22, + double freq23, + double amp, double *spec, int m0, int m1, diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index 7b076559c..2f495a114 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -181,7 +181,7 @@ def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] @cython.boundscheck(False) @cython.wraparound(False) def triangle_interpolation1D(vector, np.ndarray[double, ndim=1] spectrum_amp, - double amp=1, str type="linear"): + double amp=1.0, str type="linear"): r"""Given a vector of three points, this method interpolates the between the points to form a triangle. The height of the triangle is given as `2.0/(f[2]-f[1])` where `f` is the array `vector` sorted in an ascending @@ -197,24 +197,22 @@ def triangle_interpolation1D(vector, np.ndarray[double, ndim=1] spectrum_amp, default value is 0. :ivar type: Linear or Gaussian interpolation for delta functions. """ - cdef np.ndarray[int, ndim=1] points = np.asarray([spectrum_amp.size/2], dtype=np.int32) + cdef int points = np.asarray([spectrum_amp.size/2], dtype=np.int32)[0] cdef np.ndarray[double, ndim=1] f_vector = np.asarray(vector, dtype=np.float64) - cdef double *f1 = &f_vector[0] - cdef double *f2 = &f_vector[1] - cdef double *f3 = &f_vector[2] - - cdef np.ndarray[double, ndim=1] amp_ = np.asarray([amp]) + cdef double f1 = f_vector[0] + cdef double f2 = f_vector[1] + cdef double f3 = f_vector[2] iso_intrp = 0 if type == "linear" else 1 - clib.triangle_interpolation1D(f1, f2, f3, &_[0], &spectrum_amp[0], - &points[0], iso_intrp) + clib.triangle_interpolation1D(f1, f2, f3, amp, &spectrum_amp[0], + points, iso_intrp) @cython.boundscheck(False) @cython.wraparound(False) def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectrum_amp, - double amp=1, str type="linear"): + double amp=1.0, str type="linear"): r"""Given a vector of three points, this method interpolates the between the points to form a triangle. The height of the triangle is given as `2.0/(f[2]-f[1])` where `f` is the array `vector` sorted in an ascending @@ -229,18 +227,16 @@ def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectr cdef np.ndarray[double, ndim=1] f1_vector = np.asarray(vector1, dtype=np.float64) cdef np.ndarray[double, ndim=1] f2_vector = np.asarray(vector2, dtype=np.float64) - cdef double *f11 = &f1_vector[0] - cdef double *f12 = &f1_vector[1] - cdef double *f13 = &f1_vector[2] - - cdef double *f21 = &f2_vector[0] - cdef double *f22 = &f2_vector[1] - cdef double *f23 = &f2_vector[2] + cdef double f11 = f1_vector[0] + cdef double f12 = f1_vector[1] + cdef double f13 = f1_vector[2] - cdef np.ndarray[double, ndim=1] amp_ = np.asarray([amp]) + cdef double f21 = f2_vector[0] + cdef double f22 = f2_vector[1] + cdef double f23 = f2_vector[2] iso_intrp = 0 if type == "linear" else 1 - clib.triangle_interpolation2D(f11, f12, f13, f21, f22, f23, &_[0], + clib.triangle_interpolation2D(f11, f12, f13, f21, f22, f23, amp, &spectrum_amp[0, 0], shape[0], shape[1], iso_intrp) @cython.boundscheck(False) diff --git a/tests/test_orientations_and_triangle_interpolation.py b/tests/test_orientations_and_triangle_interpolation.py index 4ea7818bf..fc93e3d65 100644 --- a/tests/test_orientations_and_triangle_interpolation.py +++ b/tests/test_orientations_and_triangle_interpolation.py @@ -200,9 +200,8 @@ def test_triangle_interpolation(): # plt.title("1D interpolation, span(0, 100)") # plt.legend() # plt.tight_layout() - # plt.show() # plt.savefig(f"figs/fig_1D_{i}_{scl}.pdf") - # plt.figure().clear() + # plt.close() assert np.allclose(amp_py, amp_c.real, atol=1e-15) From a824c11ad089cd1f98e53f7d11877beb0d038ad8 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Wed, 28 Aug 2024 04:41:50 -0400 Subject: [PATCH 20/27] check c type consistency --- src/c_lib/include/interpolation.h | 18 ++++++-------- src/c_lib/lib/interpolation.c | 41 +++++++++++++++---------------- src/c_lib/test/test.pxd | 2 +- src/c_lib/test/test.pyx | 2 +- 4 files changed, 30 insertions(+), 33 deletions(-) diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index bde407a74..91a91f992 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -42,23 +42,21 @@ extern void two_d_averaging(double *spec, const unsigned int freq_size, double * bool user_defined, bool interpolation); extern void hist1d(double *spec, const unsigned int freq_size, double *freq, - double *amp, int m, const unsigned int nt); + double *amp, int m); extern void hist2d(double *spec, const unsigned int freq_size, double *freq_1, - double *freq_2, double *amp, int amp_stride, int m0, int m1, - const unsigned int nt); + double *freq_2, double *amp, int amp_stride, int m0, int m1); extern void generic_2d_triangle_average(double *spec, const unsigned int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, int m0, int m1, const unsigned int position_size, - int32_t *positions, const unsigned int nt, - unsigned int iso_intrp); + int32_t *positions, unsigned int iso_intrp); extern void generic_1d_triangle_average(double *spec, const unsigned int freq_size, double *freq, double *amp, int m, const unsigned int position_size, - int32_t *positions, const unsigned int nt); + int32_t *positions); extern void triangle_interpolation1D_gaussian(double f1, double f2, double f3, double amp, double *spec, int m0); @@ -95,12 +93,12 @@ extern void triangle_interpolation2D(double f11, double f12, double f13, double * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, - int stride, int n_spec, double *spec, + const unsigned int stride, int n_spec, double *spec, unsigned int iso_intrp); extern void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, int stride, int m); + double *amp, const unsigned int stride, int m); extern void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - int nt, double *amp, int stride, int m0, int m1, - unsigned int iso_intrp); + const unsigned int nt, double *amp, int stride, + int m0, int m1, unsigned int iso_intrp); diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index ac807bbe8..998be404e 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -814,9 +814,9 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // return [clipped_line] void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, - int stride, int n_spec, double *spec, + const unsigned int stride, int n_spec, double *spec, unsigned int iso_intrp) { - int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; + unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; unsigned int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address; @@ -845,8 +845,8 @@ void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *am } void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, int stride, int m) { - int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; + double *amp, const unsigned int stride, int m) { + unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; unsigned int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq_address; @@ -898,7 +898,7 @@ void generic_1d_triangle_interpolation(double *spec, const unsigned int freq_siz } void hist1d(double *spec, const unsigned int freq_size, double *freq, double *amp, - int m, const unsigned int nt) { + int m) { unsigned int i = 0, ix; double temp_freq; @@ -920,23 +920,22 @@ void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, octahedronInterpolation(spec, freq, nt, amp_real, 1, dimension_count); octahedronInterpolation(spec + 1, freq, nt, amp_imag, 1, dimension_count); } else { - hist1d(spec, freq_size, freq, amp_real, dimension_count, nt); - hist1d(spec + 1, freq_size, freq, amp_imag, dimension_count, nt); + hist1d(spec, freq_size, freq, amp_real, dimension_count); + hist1d(spec + 1, freq_size, freq, amp_imag, dimension_count); } } else { generic_1d_triangle_average(spec, freq_size, freq, amp_real, dimension_count, - position_size, positions, nt); + position_size, positions); generic_1d_triangle_average(spec + 1, freq_size, freq, amp_imag, dimension_count, - position_size, positions, nt); + position_size, positions); } } void generic_1d_triangle_average(double *spec, const unsigned int freq_size, double *freq, double *amp, int m, - const unsigned int position_size, int32_t *positions, - const unsigned int nt) { + const unsigned int position_size, int32_t *positions) { if (positions == NULL) { - hist1d(spec, freq_size, freq, amp, m, nt); + hist1d(spec, freq_size, freq, amp, m); } else { generic_1d_triangle_interpolation(spec, freq_size, freq, amp, m, position_size, positions); @@ -955,12 +954,12 @@ void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, dimension0_count, dimension1_count, iso_intrp); } else { hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, dimension0_count, - dimension1_count, nt); + dimension1_count); } } else { generic_2d_triangle_average(spec, freq_size, freq1, freq2, amp, amp_stride, dimension0_count, dimension1_count, position_size, - positions, nt, iso_intrp); + positions, iso_intrp); } } @@ -985,7 +984,7 @@ void generic_2d_triangle_interpolation(double *spec, const unsigned int freq_siz } void hist2d(double *spec, const unsigned int freq_size, double *freq1, double *freq2, - double *amp, int amp_stride, int m0, int m1, const unsigned int nt) { + double *amp, int amp_stride, int m0, int m1) { unsigned int i = 0, ix, iy, hist_index; double temp_freq_1, temp_freq_2; @@ -1005,19 +1004,19 @@ void generic_2d_triangle_average(double *spec, const unsigned int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, int m0, int m1, const unsigned int position_size, int32_t *positions, - const unsigned int nt, unsigned int iso_intrp) { + unsigned int iso_intrp) { if (positions == NULL) { - hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1, nt); + hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1); } else { generic_2d_triangle_interpolation(spec, freq_size, freq1, freq2, amp, amp_stride, position_size, positions, m0, m1, iso_intrp); } } -void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, int nt, - double *amp, int stride, int m0, int m1, - unsigned int iso_intrp) { - int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; +void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, + const unsigned int nt, double *amp, int stride, int m0, + int m1, unsigned int iso_intrp) { + unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; unsigned int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq1_address, *freq2_address; diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index c7bda5af7..8a17725ee 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -126,7 +126,7 @@ cdef extern from "interpolation.h": double *freq, int nt, double *amp, - int stride, + const unsigned int stride, int m) cdef extern from "mrsimulator.h": diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index 2f495a114..40bce4126 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -171,7 +171,7 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): @cython.boundscheck(False) @cython.wraparound(False) -def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, int stride=1): +def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, const unsigned int stride=1): cdef int i cdef int number_of_sidebands = int(amp.shape[0]) for i in range(number_of_sidebands): From 0c767884bc4cc5a9542f819a3e4fad65d065f4a9 Mon Sep 17 00:00:00 2001 From: deepanshs Date: Wed, 28 Aug 2024 05:26:24 -0400 Subject: [PATCH 21/27] cleanup --- setup.py | 2 ++ src/c_lib/base/base_model.pxd | 8 ++++---- src/c_lib/include/interpolation.h | 5 +++-- src/c_lib/include/mrsimulator.h | 3 +-- src/c_lib/include/simulation.h | 8 ++++---- src/c_lib/lib/interpolation.c | 21 ++++++++++----------- src/c_lib/lib/method.c | 5 ++--- src/c_lib/lib/mrsimulator.c | 3 +-- src/c_lib/lib/simulation.c | 8 ++++---- src/c_lib/sandbox/sandbox.pxd | 2 +- src/c_lib/sandbox/sandbox.pyx | 2 +- src/c_lib/test/test.pxd | 2 +- 12 files changed, 34 insertions(+), 35 deletions(-) diff --git a/setup.py b/setup.py index 2e8c8d363..5bdb853a7 100644 --- a/setup.py +++ b/setup.py @@ -148,6 +148,8 @@ def __init__(self): "-O3", "-ffast-math", "-fcommon", + "-Wall", + "-Wextra", # "-msse4.2", # "-ftree-vectorize", # "-fopt-info-vec-all", diff --git a/src/c_lib/base/base_model.pxd b/src/c_lib/base/base_model.pxd index 9f8e14a48..e42cfa61a 100644 --- a/src/c_lib/base/base_model.pxd +++ b/src/c_lib/base/base_model.pxd @@ -63,7 +63,7 @@ cdef extern from "mrsimulator.h": MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, double rotor_frequency_in_Hz, - double rotor_angle_in_rad, double increment, + double rotor_angle_in_rad, bool_t allow_4th_rank) void MRS_free_plan(MRS_plan *plan) void MRS_get_amplitudes_from_plan(MRS_plan *plan, bool_t refresh) @@ -133,9 +133,9 @@ cdef extern from "simulation.h": void mrsimulator_core( # spectrum information and related amplitude double *spec, - double spectral_start, - double spectral_increment, - int number_of_points, + # double spectral_start, + # double spectral_increment, + # int number_of_points, site_struct *sites, coupling_struct *couplings, diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index 91a91f992..450dbf8e3 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -100,5 +100,6 @@ extern void octahedronInterpolation(double *spec, double *freq, const unsigned i double *amp, const unsigned int stride, int m); extern void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - const unsigned int nt, double *amp, int stride, - int m0, int m1, unsigned int iso_intrp); + const unsigned int nt, double *amp, + const unsigned int stride, int m0, int m1, + unsigned int iso_intrp); diff --git a/src/c_lib/include/mrsimulator.h b/src/c_lib/include/mrsimulator.h index 4f536def2..2a57b2c56 100644 --- a/src/c_lib/include/mrsimulator.h +++ b/src/c_lib/include/mrsimulator.h @@ -83,7 +83,6 @@ typedef struct MRS_plan MRS_plan; * @param rotor_frequency_in_Hz The sample rotation frequency in Hz. * @param rotor_angle_in_rad The polar angle in radians with respect to the * z-axis describing the axis of rotation. - * @param increment The increment along the spectroscopic dimension in Hz. * @param allow_4th_rank When true, the plan calculates matrices for * processing the fourth-rank tensors. * @return A pointer to the MRS_plan. @@ -91,7 +90,7 @@ typedef struct MRS_plan MRS_plan; MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, - double increment, bool allow_4th_rank); + bool allow_4th_rank); /** * @brief Release the memory allocated for the given mrsimulator plan. diff --git a/src/c_lib/include/simulation.h b/src/c_lib/include/simulation.h index 742fcfa4c..f2a14f909 100644 --- a/src/c_lib/include/simulation.h +++ b/src/c_lib/include/simulation.h @@ -14,10 +14,10 @@ extern void mrsimulator_core( // spectrum information and related amplitude - double *spec, // Pointer to the spectrum array (complex). - double spectral_start, // The start of the frequency spectrum. - double spectral_increment, // The increment of the frequency spectrum. - int number_of_points, // Number of points on the frequency spectrum. + double *spec, // Pointer to the spectrum array (complex). + // double spectral_start, // The start of the frequency spectrum. + // double spectral_increment, // The increment of the frequency spectrum. + // int number_of_points, // Number of points on the frequency spectrum. site_struct *sites, // Pointer to sites within a spin system. coupling_struct *couplings, // Pointer to couplings within spin system. MRS_dimension *dimensions, // Pointer to the spectral dimension. diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index 998be404e..e4c82600c 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -735,7 +735,7 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // for sec in self.clip(0, 1, 1, 0): // v1 = Point(sec[0], sec[1]) // v0 = Point(sec[2], sec[3]) -// if abs(v0.x - 1) < eps and abs(v1.x - 1) < eps \ +// if abs(v0.x - 1) < eps and abs(v1.x - 1) < eps // or abs(v0.y - 1) < eps and abs(v1.y - 1) < eps: // continue @@ -879,8 +879,7 @@ void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, } } -void generic_1d_triangle_interpolation(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, +void generic_1d_triangle_interpolation(double *spec, double *freq, double *amp, int m, const unsigned int position_size, int32_t *positions) { unsigned int pos_size = position_size; @@ -937,8 +936,7 @@ void generic_1d_triangle_average(double *spec, const unsigned int freq_size, if (positions == NULL) { hist1d(spec, freq_size, freq, amp, m); } else { - generic_1d_triangle_interpolation(spec, freq_size, freq, amp, m, position_size, - positions); + generic_1d_triangle_interpolation(spec, freq, amp, m, position_size, positions); } } @@ -963,9 +961,9 @@ void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, } } -void generic_2d_triangle_interpolation(double *spec, const unsigned int freq_size, - double *freq1, double *freq2, double *amp, - int amp_stride, const unsigned int position_size, +void generic_2d_triangle_interpolation(double *spec, double *freq1, double *freq2, + double *amp, int amp_stride, + const unsigned int position_size, int32_t *positions, int m0, int m1, unsigned int iso_intrp) { unsigned int pos_size = position_size; @@ -1008,14 +1006,15 @@ void generic_2d_triangle_average(double *spec, const unsigned int freq_size, if (positions == NULL) { hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1); } else { - generic_2d_triangle_interpolation(spec, freq_size, freq1, freq2, amp, amp_stride, + generic_2d_triangle_interpolation(spec, freq1, freq2, amp, amp_stride, position_size, positions, m0, m1, iso_intrp); } } void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - const unsigned int nt, double *amp, int stride, int m0, - int m1, unsigned int iso_intrp) { + const unsigned int nt, double *amp, + const unsigned int stride, int m0, int m1, + unsigned int iso_intrp) { unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; unsigned int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq1_address, *freq2_address; diff --git a/src/c_lib/lib/method.c b/src/c_lib/lib/method.c index f521bf3d4..477e07d95 100644 --- a/src/c_lib/lib/method.c +++ b/src/c_lib/lib/method.c @@ -155,9 +155,8 @@ static inline void create_plans_for_events_in_dimension( dim->normalize_offset = 0.5 - (coordinates_offset * dim->inverse_increment); dim->R0_offset = 0.0; - MRS_plan *plan = - MRS_create_plan(scheme, number_of_sidebands, *rotor_frequency_in_Hz, - *rotor_angle_in_rad, increment, scheme->allow_4th_rank); + MRS_plan *plan = MRS_create_plan(scheme, number_of_sidebands, *rotor_frequency_in_Hz, + *rotor_angle_in_rad, scheme->allow_4th_rank); // normalize the sideband frequencies. cblas_dscal(number_of_sidebands, dim->inverse_increment, plan->vr_freq, 1); diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index 29ffefc8f..707e4ed35 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -84,7 +84,7 @@ void MRS_plan_release_temp_storage(MRS_plan *the_plan) { MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, - double increment, bool allow_4th_rank) { + bool allow_4th_rank) { MRS_plan *plan = malloc(sizeof(MRS_plan)); plan->number_of_sidebands = number_of_sidebands; plan->rotor_frequency_in_Hz = rotor_frequency_in_Hz; @@ -136,7 +136,6 @@ MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, double rotor_frequency_in_Hz) { unsigned int size_4; - // double increment_inverse = 1.0 / increment; plan->rotor_frequency_in_Hz = rotor_frequency_in_Hz; plan->is_static = (rotor_frequency_in_Hz < 1.0e-3) ? true : false; diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index 11f96e6aa..323926214 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -181,10 +181,10 @@ void __mrsimulator_core( void mrsimulator_core( // spectrum information and related amplitude - double *spec, // Pointer to the spectrum array (complex). - double coordinates_offset, // The start of the frequency spectrum. - double increment, // The increment of the frequency spectrum. - int count, // Number of points on the frequency spectrum. + double *spec, // Pointer to the spectrum array (complex). + // double coordinates_offset, // The start of the frequency spectrum. + // double increment, // The increment of the frequency spectrum. + // int count, // Number of points on the frequency spectrum. site_struct *sites, // Pointer to a list of sites within a spin system. coupling_struct *couplings, // Pointer to a list of couplings within a spin system. MRS_dimension *dimensions, // the dimensions in the method. diff --git a/src/c_lib/sandbox/sandbox.pxd b/src/c_lib/sandbox/sandbox.pxd index ff2d672de..09dd85ee1 100644 --- a/src/c_lib/sandbox/sandbox.pxd +++ b/src/c_lib/sandbox/sandbox.pxd @@ -42,7 +42,7 @@ cdef extern from "mrsimulator.h": MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, double rotor_frequency_in_Hz, - double rotor_angle_in_rad, double increment, + double rotor_angle_in_rad, bool_t allow_4th_rank) void MRS_free_plan(MRS_plan *plan) void MRS_get_amplitudes_from_plan(MRS_plan *plan, double complex *R2, diff --git a/src/c_lib/sandbox/sandbox.pyx b/src/c_lib/sandbox/sandbox.pyx index cf0d9efce..519c3b0d1 100644 --- a/src/c_lib/sandbox/sandbox.pyx +++ b/src/c_lib/sandbox/sandbox.pyx @@ -168,7 +168,7 @@ cdef class MRSPlan: self.plan = clib.MRS_create_plan(averaging_scheme.scheme, number_of_sidebands, rotor_frequency_in_Hz, - rotor_angle_in_rad, increment, allow_4th_rank) + rotor_angle_in_rad, allow_4th_rank) @property def number_of_sidebands(self): diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index 8a17725ee..e9281938f 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -141,7 +141,7 @@ cdef extern from "mrsimulator.h": # unsigned int integration_density, # int number_of_sidebands, # double rotor_frequency_in_Hz, -# double rotor_angle_in_rad, double increment, +# double rotor_angle_in_rad, # bool_t allow_4th_rank) From 6657592c40cb0997998b00ae2f24836d35ffbf70 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Wed, 28 Aug 2024 05:55:36 -0400 Subject: [PATCH 22/27] update --- setup.py | 4 +- src/c_lib/include/vm.h | 6 +-- src/mrsimulator/simulator/__init__.py | 69 ++++++++++++++------------- 3 files changed, 41 insertions(+), 38 deletions(-) diff --git a/setup.py b/setup.py index 5bdb853a7..f753be53b 100644 --- a/setup.py +++ b/setup.py @@ -185,8 +185,10 @@ def __init__(self): # "-Rpass=loop-vectorize", # "-Rpass-missed=loop-vectorize", # "-Rpass-analysis=loop-vectorize", - # "-fvectorize", + "-fvectorize", "-fcommon", + "-Wall", + "-Wextra", ] self.extra_link_args += ["-lm"] diff --git a/src/c_lib/include/vm.h b/src/c_lib/include/vm.h index b4156a17a..6d2500501 100644 --- a/src/c_lib/include/vm.h +++ b/src/c_lib/include/vm.h @@ -18,9 +18,9 @@ // absolute double value static inline double absd(double a) { - // *((unsigned __int64_ *)&a) &= ~(1ULL << 63); - // return a; - return fabs(a); + *((unsigned __int64_ *)&a) &= ~(1ULL << 63); + return a; + // return fabs(a); } /** Arithmetic suit ======================================================== */ diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 82ff585d2..aad29de34 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -8,6 +8,8 @@ import numpy as np import pandas as pd import psutil +from joblib import delayed +from joblib import Parallel from mrsimulator import Site from mrsimulator import SpinSystem from mrsimulator.base_model import core_simulator @@ -20,9 +22,6 @@ from .config import ConfigSimulator -# from joblib import delayed -# from joblib import Parallel - __author__ = "Deepansh Srivastava" __email__ = "srivastava.89@osu.edu" @@ -384,7 +383,7 @@ def run( >>> sim.run() # doctest:+SKIP """ - # verbose = 0 + verbose = 0 if opt is None: opt = self.optimize() @@ -398,37 +397,39 @@ def run( method = self.methods[index] config_dict = self.config.get_int_dict() - amp = core_simulator( - method=method, - spin_systems=self.spin_systems, - transition_pathways=opt["precomputed_pathways"][index], - transition_weights=opt["precomputed_weights"][index], - **config_dict, - **kwargs, - ) - - # spin_sys = get_chunks(self.spin_systems, n_jobs) - - # # Chunk transition pathways and weights form the optimization dictionary - # pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) - # weights = get_chunks(opt["precomputed_weights"][index], n_jobs) - - # jobs = ( - # delayed(core_simulator)( - # method=method, - # spin_systems=sys, - # transition_pathways=pth, - # transition_weights=wht, - # **config_dict, - # **kwargs, - # ) - # for sys, pth, wht in zip(spin_sys, pathways, weights) + # amp = core_simulator( + # method=method, + # spin_systems=self.spin_systems, + # transition_pathways=opt["precomputed_pathways"][index], + # transition_weights=opt["precomputed_weights"][index], + # **config_dict, + # **kwargs, # ) - # amp = Parallel( - # n_jobs=n_jobs, - # verbose=verbose, - # backend="loky", - # )(jobs) + + # Chunk spin systems for multiple jobs + spin_sys = get_chunks(self.spin_systems, n_jobs) + + # Chunk transition pathways and weights form the optimization dictionary + pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) + weights = get_chunks(opt["precomputed_weights"][index], n_jobs) + + jobs = ( + delayed(core_simulator)( + method=method, + spin_systems=sys, + transition_pathways=pth, + transition_weights=wht, + **config_dict, + **kwargs, + ) + for sys, pth, wht in zip(spin_sys, pathways, weights) + ) + amp = Parallel( + n_jobs=n_jobs, + verbose=verbose, + backend="loky", + )(jobs) + amp = amp[0] gyromagnetic_ratio = method.channels[0].gyromagnetic_ratio B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density From fc0758601db6c3ecadfe6432a235e5a224b99370 Mon Sep 17 00:00:00 2001 From: deepanshs Date: Wed, 28 Aug 2024 19:11:47 -0400 Subject: [PATCH 23/27] change all uint to int --- setup.py | 3 + src/c_lib/base/base_model.pxd | 36 ++++----- src/c_lib/base/base_model.pyx | 32 ++++---- .../include/angular_momentum/wigner_matrix.h | 7 +- src/c_lib/include/frequency_averaging.h | 7 +- src/c_lib/include/interpolation.h | 53 ++++++------- src/c_lib/include/method.h | 7 +- src/c_lib/include/mrsimulator.h | 13 ++-- src/c_lib/include/object_struct.h | 4 +- src/c_lib/include/octahedron.h | 14 ++-- src/c_lib/include/schemes.h | 57 +++++++------- src/c_lib/include/simulation.h | 12 +-- src/c_lib/include/vm_common.h | 8 +- .../lib/angular_momentum/wigner_element.c | 18 +++-- .../lib/angular_momentum/wigner_matrix.c | 22 +++--- src/c_lib/lib/method.c | 30 ++++---- src/c_lib/lib/octahedron.c | 43 +++++------ src/c_lib/lib/schemes.c | 76 +++++++++---------- src/c_lib/sandbox/sandbox.pxd | 20 ++--- src/c_lib/test/test.pxd | 24 +++--- src/c_lib/test/test.pyx | 16 ++-- 21 files changed, 248 insertions(+), 254 deletions(-) diff --git a/setup.py b/setup.py index f753be53b..a1734898f 100644 --- a/setup.py +++ b/setup.py @@ -150,6 +150,9 @@ def __init__(self): "-fcommon", "-Wall", "-Wextra", + "-Wconversion", + # "-Werror", + # "-pedantic", # "-msse4.2", # "-ftree-vectorize", # "-fopt-info-vec-all", diff --git a/src/c_lib/base/base_model.pxd b/src/c_lib/base/base_model.pxd index e42cfa61a..5b14cf4a0 100644 --- a/src/c_lib/base/base_model.pxd +++ b/src/c_lib/base/base_model.pxd @@ -26,42 +26,42 @@ cdef extern from "angular_momentum/wigner_matrix.h": cdef extern from "schemes.h": ctypedef struct MRS_averaging_scheme: - unsigned int total_orientations + int total_orientations ctypedef struct MRS_fftw_scheme: pass MRS_averaging_scheme *MRS_create_averaging_scheme( - unsigned int integration_density, + int integration_density, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, + int n_gamma, + int integration_volume, bool_t interpolation) MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( double *alpha, double *beta, double *weight, - unsigned int n_angles, + int n_angles, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int position_size, + int n_gamma, + int position_size, int32_t *positions, bool_t interpolation) void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) - MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands) + MRS_fftw_scheme *create_fftw_scheme(int total_orientations, + int number_of_sidebands) void MRS_free_fftw_scheme(MRS_fftw_scheme *fftw_scheme) cdef extern from "mrsimulator.h": ctypedef struct MRS_plan: - unsigned int number_of_sidebands + int number_of_sidebands double rotor_frequency_in_Hz double rotor_angle_in_rad - MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, + MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, bool_t allow_4th_rank) @@ -109,7 +109,7 @@ cdef extern from "method.h": double increment # Increment of coordinates along the dimension. double coordinates_offset # Start coordinate of the dimension. MRS_event *events # Holds a list of events. - unsigned int n_events # The number of events. + int n_events # The number of events. MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, @@ -123,8 +123,8 @@ cdef extern from "method.h": double *rotor_frequency_in_Hz, double *rotor_angle_in_rad, int *n_events, - unsigned int n_dim, - unsigned int *number_of_sidebands) + int n_dim, + int *number_of_sidebands) void MRS_free_dimension(MRS_dimension *dimensions, int n) @@ -146,14 +146,14 @@ cdef extern from "simulation.h": int quad_second_order, # Quad theory for second order, # spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, + int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, float *transition_pathway, # Pointer to a list of transitions. int integration_density, - unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere - unsigned int interpolate_type, + int integration_volume, # 0-octant, 1-hemisphere, 2-sphere + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix, ) @@ -169,7 +169,7 @@ cdef extern from "simulation.h": MRS_dimension *dimensions, # the dimensions within method. MRS_fftw_scheme *fftw_scheme, # the fftw scheme MRS_averaging_scheme *scheme, # the powder averaging scheme - unsigned int interpolate_type, + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix, ) diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index f90944c94..7fed0ff87 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -19,12 +19,12 @@ def core_simulator(method, list transition_pathways, # Same length as spin_systems list and list transition_weights, # elements coorespond to each system int verbose=0, # for debug purpose only. - unsigned int number_of_sidebands=90, - unsigned int integration_density=72, - unsigned int decompose_spectrum=0, - unsigned int integration_volume=1, - unsigned int isotropic_interpolation=0, - unsigned int number_of_gamma_angles=1, + int number_of_sidebands=90, + int integration_density=72, + int decompose_spectrum=0, + int integration_volume=1, + int isotropic_interpolation=0, + int number_of_gamma_angles=1, bool_t interpolation=True, bool_t auto_switch=True, debug=False, @@ -56,15 +56,15 @@ def core_simulator(method, # create averaging scheme _____________________________________________________ cdef clib.MRS_averaging_scheme *averaging_scheme - cdef unsigned int position_size = 0 + cdef int position_size = 0 if user_defined: if interpolation: - position_size = np.uint32(positions.size / 3) if positions is not None else 0 + position_size = np.intc(positions.size / 3) if positions is not None else 0 else: positions = None averaging_scheme = clib.MRS_create_averaging_scheme_from_alpha_beta( - alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=np.uint32(alpha.size), + alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=np.intc(alpha.size), allow_4th_rank=allow_4th_rank, n_gamma=number_of_gamma_angles, position_size=position_size, positions=&positions[0], interpolation=interpolation ) @@ -91,7 +91,7 @@ def core_simulator(method, cdef ndarray[int] cnt cdef ndarray[double] coord_off cdef ndarray[double] incre - cdef ndarray[unsigned int] n_dim_sidebands + cdef ndarray[int] n_dim_sidebands # Loop through dimensions and grab attributes/values in python freq_contrib = np.asarray([]) @@ -155,11 +155,11 @@ def core_simulator(method, magnetic_flux_density_in_T = np.asarray(Bo, dtype=np.float64) srfiH = np.asarray(vr, dtype=np.float64) rair = np.asarray(th, dtype=np.float64) - cnt = np.asarray(count, dtype=np.int32) + cnt = np.asarray(count, dtype=np.intc) incre = np.asarray(increment, dtype=np.float64) coord_off = np.asarray(coordinates_offset, dtype=np.float64) - n_event = np.asarray(event_i, dtype=np.int32) - n_dim_sidebands = np.asarray(dim_sidebands, dtype=np.uint32) + n_event = np.asarray(event_i, dtype=np.intc) + n_dim_sidebands = np.asarray(dim_sidebands, dtype=np.intc) # # special 1D case with 1 event. # if np.all(srfiH == 1e-3) and np.all(rair - rair[0] == 0): @@ -182,7 +182,7 @@ def core_simulator(method, norm = np.abs(np.prod(incre)) # create fftw scheme __________________________________________________________ - cdef unsigned int max_sidebands = n_dim_sidebands.max() + cdef int max_sidebands = n_dim_sidebands.max() cdef clib.MRS_fftw_scheme *fftw_scheme fftw_scheme = clib.create_fftw_scheme(averaging_scheme.total_orientations, max_sidebands) # _____________________________________________________________________________ @@ -206,7 +206,7 @@ def core_simulator(method, affine_matrix_c[3] -= affine_matrix_c[1]*affine_matrix_c[2] # sites _______________________________________________________________________________ - cdef unsigned int number_of_sites, number_of_couplings, _site_idx_, _coup_idx_ + cdef int number_of_sites, number_of_couplings, _site_idx_, _coup_idx_ cdef ndarray[int] spin_index_ij cdef ndarray[float] spin_i cdef ndarray[double] gyromagnetic_ratio_i @@ -346,7 +346,7 @@ def core_simulator(method, couplings_c.number_of_couplings = 0 if spin_sys.couplings is not None: number_of_couplings = int(len(spin_sys.couplings)) - spin_index_ij = np.zeros(2*number_of_couplings, dtype=np.int32) + spin_index_ij = np.zeros(2*number_of_couplings, dtype=np.intc) iso_j = np.zeros(number_of_couplings, dtype=np.float64) zeta_j = np.zeros(number_of_couplings, dtype=np.float64) diff --git a/src/c_lib/include/angular_momentum/wigner_matrix.h b/src/c_lib/include/angular_momentum/wigner_matrix.h index 0372d52c2..4d4234c72 100644 --- a/src/c_lib/include/angular_momentum/wigner_matrix.h +++ b/src/c_lib/include/angular_momentum/wigner_matrix.h @@ -121,8 +121,7 @@ extern void single_wigner_rotation(const int l, const double *euler_angles, * fourth-rank wigner matrices. The length of w4 is `octant_orientations x * n_octants x 9` with 9 as the leading dimension. */ -extern void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, +extern void __batch_wigner_rotation(const int octant_orientations, const int n_octants, double *wigner_2j_matrices, complex128 *R2, double *wigner_4j_matrices, complex128 *R4, complex128 *exp_Im_alpha, complex128 *w2, @@ -133,5 +132,5 @@ extern void __batch_wigner_rotation(const unsigned int octant_orientations, * The function accepts cos_angle = cos(angle). * The result is stored in exp_Im_angle as m x n matrix where m = [-4,-3,-2,-1] */ -extern void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, - void *exp_Im_angle, double delta_alpha); +extern void get_exp_Im_angle(const int n, const bool allow_4th_rank, void *exp_Im_angle, + double delta_alpha); diff --git a/src/c_lib/include/frequency_averaging.h b/src/c_lib/include/frequency_averaging.h index 61b7bc335..7bd85deb7 100644 --- a/src/c_lib/include/frequency_averaging.h +++ b/src/c_lib/include/frequency_averaging.h @@ -10,9 +10,8 @@ #include "simulation.h" void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, unsigned int iso_intrp, - complex128 *exp_I_phase); + double *spec, int iso_intrp, complex128 *exp_I_phase); void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, double *affine_matrix, - unsigned int iso_intrp, complex128 *exp_I_phase); + double *spec, double *affine_matrix, int iso_intrp, + complex128 *exp_I_phase); diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index 450dbf8e3..96fd86891 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -23,39 +23,35 @@ * 1. gaussian-interpolation. */ extern void triangle_interpolation1D(double f1, double f2, double f3, double amp, - double *spec, int m0, unsigned int iso_intrp); + double *spec, int m0, int iso_intrp); extern void triangle_interpolation1D_linear(double f1, double f2, double f3, double amp, double *spec, int m0); -extern void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, +extern void one_d_averaging(double *spec, const int freq_size, double *freq, double *amp_real, double *amp_imag, int dimension_count, - const unsigned int position_size, int32_t *positions, - const unsigned int nt, bool user_defined, - bool interpolation); + const int position_size, int32_t *positions, const int nt, + bool user_defined, bool interpolation); -extern void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, +extern void two_d_averaging(double *spec, const int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, - const unsigned int position_size, int32_t *positions, - int dimension0_count, int dimension1_count, - unsigned int iso_intrp, const unsigned int nt, - bool user_defined, bool interpolation); + const int position_size, int32_t *positions, + int dimension0_count, int dimension1_count, int iso_intrp, + const int nt, bool user_defined, bool interpolation); -extern void hist1d(double *spec, const unsigned int freq_size, double *freq, - double *amp, int m); +extern void hist1d(double *spec, const int freq_size, double *freq, double *amp, int m); -extern void hist2d(double *spec, const unsigned int freq_size, double *freq_1, - double *freq_2, double *amp, int amp_stride, int m0, int m1); +extern void hist2d(double *spec, const int freq_size, double *freq_1, double *freq_2, + double *amp, int amp_stride, int m0, int m1); -extern void generic_2d_triangle_average(double *spec, const unsigned int freq_size, +extern void generic_2d_triangle_average(double *spec, const int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, int m0, int m1, - const unsigned int position_size, - int32_t *positions, unsigned int iso_intrp); + const int position_size, int32_t *positions, + int iso_intrp); -extern void generic_1d_triangle_average(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, - const unsigned int position_size, +extern void generic_1d_triangle_average(double *spec, const int freq_size, double *freq, + double *amp, int m, const int position_size, int32_t *positions); extern void triangle_interpolation1D_gaussian(double f1, double f2, double f3, @@ -78,7 +74,7 @@ extern void triangle_interpolation1D_gaussian(double f1, double f2, double f3, */ extern void triangle_interpolation2D(double f11, double f12, double f13, double f21, double f22, double f23, double amp, double *spec, - int m0, int m1, unsigned int iso_intrp); + int m0, int m1, int iso_intrp); /** * @brief Sum amplitudes from the triangles interpolations over the region of an octant. @@ -92,14 +88,13 @@ extern void triangle_interpolation2D(double f11, double f12, double f13, double * @param spec A pointer to the starting index of a one-dimensional array * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ -void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, - const unsigned int stride, int n_spec, double *spec, - unsigned int iso_intrp); +void octahedronDeltaInterpolation(const int nt, double freq, double *amp, + const int stride, int n_spec, double *spec, + int iso_intrp); -extern void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, const unsigned int stride, int m); +extern void octahedronInterpolation(double *spec, double *freq, const int nt, + double *amp, const int stride, int m); extern void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - const unsigned int nt, double *amp, - const unsigned int stride, int m0, int m1, - unsigned int iso_intrp); + const int nt, double *amp, const int stride, + int m0, int m1, int iso_intrp); diff --git a/src/c_lib/include/method.h b/src/c_lib/include/method.h index 17f82ef79..862d57cf8 100644 --- a/src/c_lib/include/method.h +++ b/src/c_lib/include/method.h @@ -29,7 +29,7 @@ typedef struct MRS_dimension { double increment; /**< Increment of coordinates along the dimension. */ double coordinates_offset; /**< Start coordinate of the dimension. */ MRS_event *events; /**< Holds a list of events. */ - unsigned int n_events; /**< The number of events. */ + int n_events; /**< The number of events. */ /* private attributes */ double R0_offset; // holds the isotropic offset. This is used in determining if or @@ -68,8 +68,7 @@ MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, int *count, double *coordinates_offset, double *increment, double *fractions, double *durations, unsigned char *is_spectral, double *magnetic_flux_density_in_T, double *rotor_frequency_in_Hz, - double *rotor_angle_in_rad, int *n_events, unsigned int n_dim, - unsigned int *number_of_sidebands); + double *rotor_angle_in_rad, int *n_events, int n_dim, int *number_of_sidebands); /** * @brief Free the memory allocation for the MRS dimensions. @@ -77,6 +76,6 @@ MRS_dimension *MRS_create_dimensions( * @param dimensions The pointer to an array of MRS_dimension structs. * @param n An interger defining the number of MRS_dimension structs in `dimensions`. */ -void MRS_free_dimension(MRS_dimension *dimensions, unsigned int n); +void MRS_free_dimension(MRS_dimension *dimensions, int n); #endif /* method_h */ diff --git a/src/c_lib/include/mrsimulator.h b/src/c_lib/include/mrsimulator.h index 2a57b2c56..286f280ac 100644 --- a/src/c_lib/include/mrsimulator.h +++ b/src/c_lib/include/mrsimulator.h @@ -34,8 +34,8 @@ */ struct MRS_plan { - unsigned int number_of_sidebands; /**< The number of sidebands to compute. */ - double rotor_frequency_in_Hz; /**< The sample rotation frequency in Hz. */ + int number_of_sidebands; /**< The number of sidebands to compute. */ + double rotor_frequency_in_Hz; /**< The sample rotation frequency in Hz. */ /** * The angle, in radians, describing the sample axis-rotation with respect to @@ -62,8 +62,8 @@ struct MRS_plan { bool copy_for_rotor_freq; // Set True if plan is copied from rotor freq update. bool allow_4th_rank; // If true, creates buffer/tables for 4th-rank tensors. bool is_static; // It true, compute static frequencies - unsigned int size; // # of angular orientations * number of sizebands. - unsigned int n_octants; // # of octants used in the orientational averaging. + int size; // number of angular orientations * number of sidebands. + int n_octants; // number of octants used in the orientational averaging. double *norm_amplitudes; // array of normalized amplitudes per orientation. double *wigner_d2m0_vector; // wigner-2j dm0 vector, n ∈ [-2, 2]. double *wigner_d4m0_vector; // wigner-4j dm0 vector, n ∈ [-4, 4]. @@ -87,8 +87,7 @@ typedef struct MRS_plan MRS_plan; * processing the fourth-rank tensors. * @return A pointer to the MRS_plan. */ -MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, - unsigned int number_of_sidebands, +MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, bool allow_4th_rank); @@ -223,7 +222,7 @@ void MRS_rotate_components_from_PAS_to_common_frame( unsigned char *freq_contrib // The pointer to freq contribs boolean. ); -extern void get_sideband_phase_components(unsigned int number_of_sidebands, +extern void get_sideband_phase_components(int number_of_sidebands, double spin_frequency, double *restrict pre_phase); diff --git a/src/c_lib/include/object_struct.h b/src/c_lib/include/object_struct.h index e4ad61855..5b248a410 100644 --- a/src/c_lib/include/object_struct.h +++ b/src/c_lib/include/object_struct.h @@ -17,7 +17,7 @@ * Site structure is a collection of site interaction parameters. **/ struct __site_struct { - unsigned int number_of_sites; /**< Number of sites */ + int number_of_sites; /**< Number of sites */ /* Pointer to an array of spin quantum numbers for each site within a spin system. */ float *spin; @@ -67,7 +67,7 @@ typedef struct __site_struct site_struct; * Coupling structure is a collection of coupled site interaction parameters. **/ struct __coupling_struct { - unsigned int number_of_couplings; /**< Number of couplings */ + int number_of_couplings; /**< Number of couplings */ /* Pointer to an array of the site indexes for each coupled pair within a spin system, * incremented in stride of 2/pair. The array size is 2*number_of_couplings. */ diff --git a/src/c_lib/include/octahedron.h b/src/c_lib/include/octahedron.h index 683a9002e..96ea3041a 100644 --- a/src/c_lib/include/octahedron.h +++ b/src/c_lib/include/octahedron.h @@ -11,20 +11,20 @@ #include "interpolation.h" extern void octahedronGetDirectionCosineSquareAndWeightsOverOctant( - const unsigned int nt, double *restrict xr, double *restrict yr, - double *restrict zr, double *restrict amp); + const int nt, double *restrict xr, double *restrict yr, double *restrict zr, + double *restrict amp); -extern void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, +extern void octahedronGetPolarAngleTrigOverOctant(const int nt, double *restrict cos_alpha, double *restrict cos_beta, double *restrict amp); -extern void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, +extern void octahedronGetComplexExpOfPolarAngleOverOctant(const int nt, void *restrict exp_I_alpha, void *restrict exp_I_beta, double *restrict amp); -void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum); +void get_total_amplitude(const int nt, double *amp, double *amp_sum); -extern void averaging_setup(unsigned int nt, void *exp_I_alpha, void *exp_I_beta, - double *amp, bool interpolation); +extern void averaging_setup(int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, + bool interpolation); diff --git a/src/c_lib/include/schemes.h b/src/c_lib/include/schemes.h index 87e72ceea..869aa5d12 100644 --- a/src/c_lib/include/schemes.h +++ b/src/c_lib/include/schemes.h @@ -36,30 +36,30 @@ * thousands of sites. */ typedef struct MRS_averaging_scheme { - unsigned int total_orientations; /**< The total number of orientations. */ + int total_orientations; /**< The total number of orientations. */ /** \privatesection */ - unsigned int integration_density; // # triangles along the edge of the octahedron. - unsigned int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. - unsigned int octant_orientations; // # unique orientations on the face of an octant. - unsigned int n_gamma; // number of gamma angles - double *amplitudes; // array of amplitude scaling per orientation. - complex128 *exp_Im_alpha; // array of e^in\alpha n=[0,4] per orientation. - complex128 *exp_Im_gamma; // array of e^im\gamma m=[0,4] per orientation. - complex128 *w2; // buffer for 2nd rank frequency calculation. - complex128 *w4; // buffer for 4nd rank frequency calculation. - double *wigner_2j_matrices; // wigner-d 2j matrix per orientation. - double *wigner_4j_matrices; // wigner-d 4j matrix per orientation. - double *scratch; // scratch memory for calculations. - bool allow_4th_rank; // If true, compute wigner matrices for wigner-d 4j. - double *amps_real; // array of real amplitudes - double *amps_imag; // array of real amplitudes + int integration_density; // number of triangles along the edge of the octahedron. + int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. + int octant_orientations; // number of unique orientations on the face of an octant. + int n_gamma; // number of gamma angles + double *amplitudes; // array of amplitude scaling per orientation. + complex128 *exp_Im_alpha; // array of e^in\alpha n=[0,4] per orientation. + complex128 *exp_Im_gamma; // array of e^im\gamma m=[0,4] per orientation. + complex128 *w2; // buffer for 2nd rank frequency calculation. + complex128 *w4; // buffer for 4nd rank frequency calculation. + double *wigner_2j_matrices; // wigner-d 2j matrix per orientation. + double *wigner_4j_matrices; // wigner-d 4j matrix per orientation. + double *scratch; // scratch memory for calculations. + bool allow_4th_rank; // If true, compute wigner matrices for wigner-d 4j. + double *amps_real; // array of real amplitudes + double *amps_imag; // array of real amplitudes double *phase; complex128 *exp_I_phase; - int32_t *positions; // array of triangle vertexes (faces) on mesh. - unsigned int position_size; // number of triangle vertexes (faces) on mesh - bool user_defined; // if true, the scheme is user defined - bool interpolation; // if true, use frequency triangle interpolation + int32_t *positions; // array of triangle vertexes (faces) on mesh. + int position_size; // number of triangle vertexes (faces) on mesh + bool user_defined; // if true, the scheme is user defined + bool interpolation; // if true, use frequency triangle interpolation } MRS_averaging_scheme; // typedef struct MRS_averaging_scheme; @@ -82,10 +82,9 @@ typedef struct MRS_averaging_scheme { * @param integration_volume An enumeration. 0=octant, 1=hemisphere, 2=sphere * @param interpolation If true, apply frequency triangle interpolation. */ -MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_density, - bool allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, +MRS_averaging_scheme *MRS_create_averaging_scheme(int integration_density, + bool allow_4th_rank, int n_gamma, + int integration_volume, bool interpolation); /** @@ -97,7 +96,7 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi * double. * @param weight A pointer to an array of size `n_angles` holding the weights of the * angle pair (alpha, beta) of type double. - * @param n_angles An unsigned int equal to the total number of alpha, beta, and weight + * @param n_angles An int equal to the total number of alpha, beta, and weight * values. * @param allow_4th_rank If true, the scheme also calculates matrices for fourth-rank * tensors. @@ -107,9 +106,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi * @param interpolation If true, apply frequency triangle interpolation. */ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( - double *alpha, double *beta, double *weight, unsigned int n_angles, - bool allow_4th_rank, unsigned int n_gamma, const unsigned int position_size, - int32_t *positions, bool interpolation); + double *alpha, double *beta, double *weight, int n_angles, bool allow_4th_rank, + int n_gamma, const int position_size, int32_t *positions, bool interpolation); /** * Free the memory allocated for the spatial orientation averaging scheme. @@ -132,8 +130,7 @@ typedef struct MRS_fftw_scheme { fftw_plan the_fftw_plan; // The plan for fftw routine. } MRS_fftw_scheme; -MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands); +MRS_fftw_scheme *create_fftw_scheme(int total_orientations, int number_of_sidebands); void MRS_free_fftw_scheme(MRS_fftw_scheme *fftw_scheme); diff --git a/src/c_lib/include/simulation.h b/src/c_lib/include/simulation.h index f2a14f909..8b8515eae 100644 --- a/src/c_lib/include/simulation.h +++ b/src/c_lib/include/simulation.h @@ -24,16 +24,16 @@ extern void mrsimulator_core( int n_dimension, // Number of dimensions. int quad_second_order, // Quad theory for second order, - unsigned int number_of_sidebands, // The number of sidebands. - double rotor_frequency_in_Hz, // The rotor spin frequency. - double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis. + int number_of_sidebands, // The number of sidebands. + double rotor_frequency_in_Hz, // The rotor spin frequency. + double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis. float *transition_pathway, // Pointer to a list of transitions. double *transition_pathway_weight, // The complex weight of transition pathway. // powder orientation average int integration_density, // The number of triangle along the edge of octahedron. - unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. - unsigned int interpolate_type, unsigned char *freq_contrib, double *affine_matrix); + int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix); extern void __mrsimulator_core( // spectrum information and related amplitude @@ -51,7 +51,7 @@ extern void __mrsimulator_core( MRS_dimension *dimensions, // Pointer to MRS_dimension structure. MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - unsigned int interpolate_type, + int interpolate_type, /** * Each event consists of the following freq contrib ordered as diff --git a/src/c_lib/include/vm_common.h b/src/c_lib/include/vm_common.h index 4a90f0383..ba7ae7c19 100644 --- a/src/c_lib/include/vm_common.h +++ b/src/c_lib/include/vm_common.h @@ -107,7 +107,7 @@ static inline void vm_double_ones(int count, double *restrict res) { * @returns values A pointer to the fft output order vector of size @p n. */ static inline double *get_FFT_order_freq(int n, double increment) { - double *vr_freq = malloc_double(n); + double *vr_freq = malloc_double((size_t)n); int m, positive_limit, negative_limit; if (n % 2 == 0) { @@ -121,7 +121,7 @@ static inline double *get_FFT_order_freq(int n, double increment) { for (m = 0; m <= positive_limit; m++) *vr_freq++ = (double)m * increment; for (m = negative_limit; m < 0; m++) *vr_freq++ = (double)m * increment; return vr_freq - n; -}; +} /** * @brief Return a vector ordered according to the fft output order. @@ -130,7 +130,7 @@ static inline double *get_FFT_order_freq(int n, double increment) { * @returns values A pointer to the fft output order vector of size @p n. */ static inline int *get_FFT_index(int n) { - int *freq_index = malloc_int(n); + int *freq_index = malloc_int((size_t)n); int m, positive_limit, negative_limit; if (n % 2 == 0) { @@ -144,4 +144,4 @@ static inline int *get_FFT_index(int n) { for (m = 0; m <= positive_limit; m++) *freq_index++ = m; for (m = negative_limit; m < 0; m++) *freq_index++ = m; return freq_index - n; -}; +} diff --git a/src/c_lib/lib/angular_momentum/wigner_element.c b/src/c_lib/lib/angular_momentum/wigner_element.c index 63a48eba5..df3e106aa 100644 --- a/src/c_lib/lib/angular_momentum/wigner_element.c +++ b/src/c_lib/lib/angular_momentum/wigner_element.c @@ -37,16 +37,17 @@ static inline double __generic_wigner_d_element(const float l, const float m1, double cx = cos(beta / 2.); double sum = 0.0, k1, k2, k3, x, y; int sign = 1; - int k; + int k, n1, n2; + float two = 2.0; for (k = 0; k <= (int)(l - m1); k++) { - k1 = (int)(l - m1 - k); - k2 = (int)(l + m2 - k); - k3 = (int)(k + m1 - m2); + k1 = (int)(l - m1 - (float)k); + k2 = (int)(l + m2 - (float)k); + k3 = (int)((float)k + m1 - m2); if (k1 >= 0 && k2 >= 0 && k3 >= 0) { - int n1 = (int)(2 * l + m2 - m1 - 2 * k); - int n2 = (int)(m1 - m2 + 2 * k); + n1 = (int)(two * l + m2 - m1 - two * (float)k); + n2 = (int)(m1 - m2 + two * (float)k); x = my_power(cx, n1); y = my_power(sx, n2); sum += sign * x * y / @@ -54,7 +55,8 @@ static inline double __generic_wigner_d_element(const float l, const float m1, } sign = -sign; } - double f = fac(l + m2) * fac(l - m2) * fac(l + m1) * fac(l - m1); + double f = fac((double)(l + m2)) * fac((double)(l - m2)) * fac((double)(l + m1)) * + fac((double)(l - m1)); sign = ((int)(m1 - m2) % 2 == 0) ? 1 : -1; f = sqrt(f); return (sign * sum * f); @@ -358,7 +360,7 @@ void transition_connect_factor(const float l, const float m1_f, const float m1_i // const double phi, double *restrict factor, int // n_sites) { // double m1_f, m1_i, m2_f, m2_i; -// complex128 *weight = malloc_complex128(n_sites); +// complex128 *weight = malloc_complex128((size_t)n_sites); // for (int i = 1; i < n_sites; i++) { // m1_i = transition_inital[i]; // m1_f = (transition_inital + n_sites)[i]; diff --git a/src/c_lib/lib/angular_momentum/wigner_matrix.c b/src/c_lib/lib/angular_momentum/wigner_matrix.c index 3f1eae6c4..3f5f044b9 100644 --- a/src/c_lib/lib/angular_momentum/wigner_matrix.c +++ b/src/c_lib/lib/angular_momentum/wigner_matrix.c @@ -17,7 +17,7 @@ complex128 NEGATIVE_IOTA = {0.0, -1.0}; // ✅ .. note: (wigner_d_matrices) tested with pytest // ......................... void wigner_d_matrices(const int l, const int n, const double *beta, double *wigner) { - complex128 *exp_I_beta = malloc_complex128(n); + complex128 *exp_I_beta = malloc_complex128((size_t)n); vm_cosine_I_sine(n, beta, exp_I_beta); wigner_d_matrices_from_exp_I_beta(l, n, false, exp_I_beta, wigner); free(exp_I_beta); @@ -741,12 +741,12 @@ void single_wigner_rotation(const int l, const double *euler_angles, const void * rotation with fourth rank wigner matrices. The length of w4 is * `octant_orientations x n_octants x 9` with 9 as the leading dimension. */ -void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, double *wigner_2j_matrices, - complex128 *R2, double *wigner_4j_matrices, complex128 *R4, +void __batch_wigner_rotation(const int octant_orientations, const int n_octants, + double *wigner_2j_matrices, complex128 *R2, + double *wigner_4j_matrices, complex128 *R4, complex128 *exp_Im_alpha, complex128 *w2, complex128 *w4) { - unsigned int j, max_iter, wigner_2j_inc, wigner_4j_inc = 0; - unsigned int w2_increment, w4_increment = 0, alpha_inc; + int j, max_iter, wigner_2j_inc, wigner_4j_inc = 0; + int w2_increment, w4_increment = 0, alpha_inc; double *exp_Im_alpha_ = (double *)exp_Im_alpha; w2_increment = 3 * octant_orientations; @@ -841,13 +841,13 @@ void __batch_wigner_rotation(const unsigned int octant_orientations, * The result is stored in exp_Im_angle as m x n matrix, where m = [-4,-3,-2,-1] * Since only negative m's are computed, the following computes exp(Im angle) */ -void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, - void *exp_Im_angle, double delta_alpha) { +void get_exp_Im_angle(const int n, const bool allow_4th_rank, void *exp_Im_angle, + double delta_alpha) { double *exp_Im_angle_ = (double *)exp_Im_angle; // The complex array is interpreted as alternating real and imag double array. // The index s_i = i*n of complex array is now at index 2*i*n. - unsigned int m_3 = 2 * n, m_2 = 4 * n, m_1 = 6 * n; // m_4 is 0 + int m_3 = 2 * n, m_2 = 4 * n, m_1 = 6 * n; // m_4 is 0 // exp(-2 I angle) vm_double_complex_multiply(n, &exp_Im_angle_[m_1], &exp_Im_angle_[m_1], @@ -863,8 +863,8 @@ void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, } if (delta_alpha != 0.0) { - double *exp_da = malloc_double(8); - unsigned int m_4p = 8 * n, m_3p = 10 * n, m_2p = 12 * n, m_1p = 14 * n; + double *exp_da = malloc_double((size_t)8); + int m_4p = 8 * n, m_3p = 10 * n, m_2p = 12 * n, m_1p = 14 * n; exp_da[0] = cos(delta_alpha); exp_da[1] = sin(delta_alpha); diff --git a/src/c_lib/lib/method.c b/src/c_lib/lib/method.c index 477e07d95..f3cd45047 100644 --- a/src/c_lib/lib/method.c +++ b/src/c_lib/lib/method.c @@ -21,8 +21,8 @@ void MRS_free_event(MRS_event *the_event) { } /** Free the memory/buffer allocation for the MRS dimensions and events within. **/ -void MRS_free_dimension(MRS_dimension *dimensions, unsigned int n) { - unsigned int dim, evt; +void MRS_free_dimension(MRS_dimension *dimensions, int n) { + int dim, evt; for (dim = 0; dim < n; dim++) { if (DEBUG) printf("post execution dimension %d cleanup\n", dim); for (evt = 0; evt < dimensions[dim].n_events; evt++) { @@ -143,13 +143,13 @@ static inline void create_plans_for_events_in_dimension( double coordinates_offset, int n_events, double *fraction, double *duration, unsigned char *is_spectral, double *rotor_frequency_in_Hz, double *rotor_angle_in_rad, double *magnetic_flux_density_in_T, - unsigned int number_of_sidebands) { + int number_of_sidebands) { int i; dim->count = count; dim->coordinates_offset = coordinates_offset; dim->increment = increment; dim->n_events = n_events; - dim->events = (MRS_event *)malloc(n_events * sizeof(MRS_event)); + dim->events = (MRS_event *)malloc((size_t)n_events * sizeof(MRS_event)); dim->inverse_increment = 1.0 / increment; dim->normalize_offset = 0.5 - (coordinates_offset * dim->inverse_increment); @@ -161,12 +161,12 @@ static inline void create_plans_for_events_in_dimension( cblas_dscal(number_of_sidebands, dim->inverse_increment, plan->vr_freq, 1); for (i = 0; i < n_events; i++) { - dim->events[i].event_freq_amplitude = malloc_complex128(plan->size); + dim->events[i].event_freq_amplitude = malloc_complex128((size_t)plan->size); vm_double_ones(plan->size * 2, (double *)dim->events[i].event_freq_amplitude); cblas_dscal(plan->size, 0.0, (double *)dim->events[i].event_freq_amplitude + 1, 2); // if (*rotor_frequency_in_Hz != 0.0 && *rotor_frequency_in_Hz != 1.0e12) { - // dim->events[i].freq_amplitude = malloc_double(plan->size); + // dim->events[i].freq_amplitude = malloc_double((size_t)plan->size); // vm_double_ones(plan->size, dim->events[i].freq_amplitude); // } MRS_set_event(&(dim->events[i]), *fraction++, *duration++, *is_spectral++, @@ -180,10 +180,12 @@ static inline void create_plans_for_events_in_dimension( /* buffer to hold the local frequencies and frequency offset. The buffer is useful * when the rotor angle is off magic angle (54.735 deg). */ - dim->local_frequency = malloc_double(scheme->n_gamma * scheme->total_orientations); - dim->local_phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - dim->freq_offset = malloc_double(scheme->octant_orientations); - dim->freq_amplitude = malloc_double(plan->size); + dim->local_frequency = + malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + dim->local_phase = + malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + dim->freq_offset = malloc_double((size_t)scheme->octant_orientations); + dim->freq_amplitude = malloc_double((size_t)plan->size); } /** @@ -194,10 +196,10 @@ MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, int *count, double *coordinates_offset, double *increment, double *fractions, double *durations, unsigned char *is_spectral, double *magnetic_flux_density_in_T, double *rotor_frequency_in_Hz, - double *rotor_angle_in_rad, int *n_events, unsigned int n_dim, - unsigned int *number_of_sidebands) { - unsigned int i; - MRS_dimension *dimension = (MRS_dimension *)malloc(n_dim * sizeof(MRS_dimension)); + double *rotor_angle_in_rad, int *n_events, int n_dim, int *number_of_sidebands) { + int i; + MRS_dimension *dimension = + (MRS_dimension *)malloc((size_t)n_dim * sizeof(MRS_dimension)); for (i = 0; i < n_dim; i++) { create_plans_for_events_in_dimension( diff --git a/src/c_lib/lib/octahedron.c b/src/c_lib/lib/octahedron.c index 3ebe0f834..72dc6e746 100644 --- a/src/c_lib/lib/octahedron.c +++ b/src/c_lib/lib/octahedron.c @@ -9,12 +9,12 @@ #include "octahedron.h" -void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const unsigned int nt, +void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const int nt, double *restrict xr, double *restrict yr, double *restrict zr, double *restrict amp) { - unsigned int i, j; + int i, j; // scale is divided by pi to remove angular dependence. // The factor 6 is because each point in the triangle interpolation is shared by // 6 neighboring triangles. @@ -48,15 +48,14 @@ void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const unsigned int n *amp = 1.0 / (r2 * r2 * factor); } -void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, - double *restrict cos_alpha, +void octahedronGetPolarAngleTrigOverOctant(const int nt, double *restrict cos_alpha, double *restrict cos_beta, double *restrict amp) { - unsigned int points = (nt + 1) * (nt + 2) / 2; - double *xr = malloc_double(points); - double *yr = malloc_double(points); - double *zr = malloc_double(points); - double *sin_beta = malloc_double(points); + int points = (nt + 1) * (nt + 2) / 2; + double *xr = malloc_double((size_t)points); + double *yr = malloc_double((size_t)points); + double *zr = malloc_double((size_t)points); + double *sin_beta = malloc_double((size_t)points); // The values xr = x^2, yr = y^2, zr = z^2, where x, y, and z are the // direction cosines. @@ -78,7 +77,7 @@ void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, // vdSqrt(points, &yr[0], &yr[0]); // Evaluate cos_alpha = x/sqrt(x^2 + y^2) = zr/sin_beta - vm_double_divide(points - 1, zr, sin_beta, cos_alpha); + vm_double_divide((points - 1), zr, sin_beta, cos_alpha); // vdDiv(points-1, yr, sinBeta, sinAlpha ); cos_alpha[points - 1] = 1.0; @@ -90,14 +89,14 @@ void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, free(sin_beta); } -void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, +void octahedronGetComplexExpOfPolarAngleOverOctant(const int nt, void *restrict exp_I_alpha, void *restrict exp_I_beta, double *restrict amp) { - unsigned int points = (nt + 1) * (nt + 2) / 2; - double *xr = malloc_double(points); - double *yr = malloc_double(points); - double *zr = malloc_double(points); + int points = (nt + 1) * (nt + 2) / 2; + double *xr = malloc_double((size_t)points); + double *yr = malloc_double((size_t)points); + double *zr = malloc_double((size_t)points); octahedronGetDirectionCosineSquareAndWeightsOverOctant(nt, xr, yr, zr, amp); // At this point the variables @@ -143,7 +142,7 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, // cos(alpha) = x/sqrt(x^2 + y^2) = x/sin(beta) ... (3) // In terms of the variables, Eq (3) is given as xr/zr. // Evaluate xr = xr/zr - vm_double_divide_inplace(points - 1, zr, xr); + vm_double_divide_inplace((points - 1), zr, xr); xr[points - 1] = 0.0; // Copy cos(alpha), aka xr, to the even addresses of exp_I_alpha cblas_dcopy(points, xr, 1, (double *)exp_I_alpha, 2); @@ -152,7 +151,7 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, // sin(alpha) = y/sqrt(x^2 + y^2) = y/sin(beta) ... (4) // In terms of the variables, Eq (4) is given as yr/zr. // Evaluate yr = yr/zr - vm_double_divide_inplace(points - 1, zr, yr); + vm_double_divide_inplace((points - 1), zr, yr); yr[points - 1] = 0.0; // Copy sin(alpha), aka yr, to the odd addresses of exp_I_alpha. cblas_dcopy(points, yr, 1, (double *)exp_I_alpha + 1, 2); @@ -163,8 +162,8 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, free(zr); } -void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum) { - unsigned int i = 0, j = 0, local_index = nt - 1, n_pts = (nt + 1) * (nt + 2) / 2; +void get_total_amplitude(const int nt, double *amp, double *amp_sum) { + int i = 0, j = 0, local_index = nt - 1, n_pts = (nt + 1) * (nt + 2) / 2; double *amp_address, temp; amp_address = &[nt + 1]; @@ -189,8 +188,8 @@ void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum) { // fix amplitude for binning case: // octant face edge points are divided by 2 // octant vertexes are divided by 4 -void fix_amplitude_for_binning(unsigned int nt, double *amp) { - unsigned int i = 0, i_max = (nt + 1) * (nt + 2) / 2, tt = nt - 1; +void fix_amplitude_for_binning(int nt, double *amp) { + int i = 0, i_max = (nt + 1) * (nt + 2) / 2, tt = nt - 1; bool inc_one = false; while (i <= nt) amp[i++] /= 2.0; while (i < i_max) { @@ -211,7 +210,7 @@ void fix_amplitude_for_binning(unsigned int nt, double *amp) { cblas_dscal(i_max, 6.0, amp, 1); } -void averaging_setup(unsigned int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, +void averaging_setup(int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, bool interpolation) { // octahedronGetPolarAngleTrigOverOctant(nt, cos_alpha, cos_beta, amp); octahedronGetComplexExpOfPolarAngleOverOctant(nt, exp_I_alpha, exp_I_beta, amp); diff --git a/src/c_lib/lib/schemes.c b/src/c_lib/lib/schemes.c index 50c343fa3..3f8485661 100644 --- a/src/c_lib/lib/schemes.c +++ b/src/c_lib/lib/schemes.c @@ -11,7 +11,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, complex128 *exp_I_beta, bool allow_4th_rank) { - unsigned int allocate_size_2, allocate_size_4; + int allocate_size_2, allocate_size_4; double delta_alpha = 0.0; scheme->total_orientations = scheme->octant_orientations; @@ -56,7 +56,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, /* Second-rank reduced wigner matrices at every β orientation from the positive upper * octant. */ - scheme->wigner_2j_matrices = malloc_double(allocate_size_2); + scheme->wigner_2j_matrices = malloc_double((size_t)allocate_size_2); wigner_d_matrices_from_exp_I_beta(2, scheme->octant_orientations, true, exp_I_beta, scheme->wigner_2j_matrices); @@ -64,7 +64,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, if (allow_4th_rank) { /* Fourt-rank reduced wigner matrices at every β orientation from the positive upper * octant. */ - scheme->wigner_4j_matrices = malloc_double(allocate_size_4); + scheme->wigner_4j_matrices = malloc_double((size_t)allocate_size_4); wigner_d_matrices_from_exp_I_beta(4, scheme->octant_orientations, true, exp_I_beta, scheme->wigner_4j_matrices); } @@ -106,13 +106,13 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, /* ................................................................................ */ /* w2 is the buffer for storing the frequencies calculated from the second-rank * tensors. Only calcuate the -2, -1, and 0 tensor components.*/ - scheme->w2 = malloc_complex128(3 * scheme->total_orientations); + scheme->w2 = malloc_complex128((size_t)(3 * scheme->total_orientations)); scheme->w4 = NULL; if (allow_4th_rank) { /* w4 is the buffer for storing the frequencies calculated from the fourth-rank * tensors. Only calcuate the -4, -3, -2, -1, and 0 tensor components.*/ - scheme->w4 = malloc_complex128(5 * scheme->total_orientations); + scheme->w4 = malloc_complex128((size_t)(5 * scheme->total_orientations)); } } @@ -135,10 +135,9 @@ void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) { } /* Create a new orientation averaging scheme. */ -MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_density, - bool allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, +MRS_averaging_scheme *MRS_create_averaging_scheme(int integration_density, + bool allow_4th_rank, int n_gamma, + int integration_volume, bool interpolation) { int alpha_size; MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); @@ -160,10 +159,10 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi // The 4 * octant_orientations memory allocation is for m=4, 3, 2, and 1 alpha_size = 4 * scheme->octant_orientations; alpha_size *= (integration_volume == 2) ? 2 : 1; - scheme->exp_Im_alpha = malloc_complex128(alpha_size); - scheme->exp_Im_gamma = malloc_complex128(4 * n_gamma); - complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); - scheme->amplitudes = malloc_double(scheme->octant_orientations); + scheme->exp_Im_alpha = malloc_complex128((size_t)alpha_size); + scheme->exp_Im_gamma = malloc_complex128((size_t)(4 * n_gamma)); + complex128 *exp_I_beta = malloc_complex128((size_t)scheme->octant_orientations); + scheme->amplitudes = malloc_double((size_t)scheme->octant_orientations); averaging_setup(integration_density, &scheme->exp_Im_alpha[3 * scheme->octant_orientations], exp_I_beta, @@ -172,8 +171,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi averaging_scheme_setup(scheme, exp_I_beta, allow_4th_rank); // exp(-im gamma) for m=[-4,-1], and gamma=[0..8]*2pi/9 - double *gamma = malloc_double(n_gamma); - double *temp = malloc_double(n_gamma); + double *gamma = malloc_double((size_t)n_gamma); + double *temp = malloc_double((size_t)n_gamma); double factor = CONST_2PI / (double)n_gamma; vm_double_arange(n_gamma, gamma); cblas_dscal(n_gamma, factor, gamma, 1); @@ -182,10 +181,11 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi // reallocate exp_I_beta memory as scratch. scheme->scratch = (double *)exp_I_beta; - scheme->amps_real = malloc_double(scheme->total_orientations); - scheme->amps_imag = malloc_double(scheme->total_orientations); - scheme->phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - scheme->exp_I_phase = malloc_complex128(scheme->n_gamma * scheme->total_orientations); + scheme->amps_real = malloc_double((size_t)scheme->total_orientations); + scheme->amps_imag = malloc_double((size_t)scheme->total_orientations); + scheme->phase = malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + scheme->exp_I_phase = + malloc_complex128((size_t)(scheme->n_gamma * scheme->total_orientations)); free(gamma); free(temp); @@ -194,9 +194,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi /* Create a new orientation averaging scheme. */ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( - double *alpha, double *beta, double *weight, unsigned int n_angles, - bool allow_4th_rank, unsigned int n_gamma, const unsigned int position_size, - int32_t *positions, bool interpolation) { + double *alpha, double *beta, double *weight, int n_angles, bool allow_4th_rank, + int n_gamma, const int position_size, int32_t *positions, bool interpolation) { double scale; MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); @@ -214,10 +213,10 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( /* Calculate α, β, and weights over the positive octant. .......................... */ /* ................................................................................ */ // The 4 * octant_orientations memory allocation is for m=4, 3, 2, and 1 - scheme->exp_Im_alpha = malloc_complex128(4 * scheme->octant_orientations); - scheme->exp_Im_gamma = malloc_complex128(4 * n_gamma); - complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); - scheme->amplitudes = malloc_double(scheme->octant_orientations); + scheme->exp_Im_alpha = malloc_complex128((size_t)(4 * scheme->octant_orientations)); + scheme->exp_Im_gamma = malloc_complex128((size_t)(4 * n_gamma)); + complex128 *exp_I_beta = malloc_complex128((size_t)scheme->octant_orientations); + scheme->amplitudes = malloc_double((size_t)scheme->octant_orientations); // scale = (1 / 6) factor for 1 point to 6 triangle ratio. scale = (interpolation) ? 0.1666666667 : 1.0; @@ -233,8 +232,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( averaging_scheme_setup(scheme, exp_I_beta, allow_4th_rank); // exp(-im gamma) for m=[-4,-1], and gamma=[0..8]*2pi/9 - double *gamma = malloc_double(n_gamma); - double *temp = malloc_double(n_gamma); + double *gamma = malloc_double((size_t)n_gamma); + double *temp = malloc_double((size_t)n_gamma); double factor = CONST_2PI / (double)n_gamma; vm_double_arange(n_gamma, gamma); cblas_dscal(n_gamma, factor, gamma, 1); @@ -243,10 +242,11 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( // reallocate exp_I_beta memory as scratch. scheme->scratch = (double *)exp_I_beta; - scheme->amps_real = malloc_double(scheme->total_orientations); - scheme->amps_imag = malloc_double(scheme->total_orientations); - scheme->phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - scheme->exp_I_phase = malloc_complex128(scheme->n_gamma * scheme->total_orientations); + scheme->amps_real = malloc_double((size_t)scheme->total_orientations); + scheme->amps_imag = malloc_double((size_t)scheme->total_orientations); + scheme->phase = malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + scheme->exp_I_phase = + malloc_complex128((size_t)(scheme->n_gamma * scheme->total_orientations)); free(gamma); free(temp); @@ -256,15 +256,15 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( /* ---------------------------------------------------------------------------------- */ /* fftw routine setup ............................................................... */ /* .................................................................................. */ -MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands) { - unsigned int size = total_orientations * number_of_sidebands; - int nssb = (int)number_of_sidebands; +MRS_fftw_scheme *create_fftw_scheme(int total_orientations, int number_of_sidebands) { + int size = total_orientations * number_of_sidebands; + int nssb = number_of_sidebands; MRS_fftw_scheme *fftw_scheme = malloc(sizeof(MRS_fftw_scheme)); // fftw_scheme->vector = fftw_alloc_complex(size); - fftw_scheme->vector = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * size); - // malloc_complex128(plan->size); + fftw_scheme->vector = + (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (size_t)size); + // malloc_complex128((size_t)plan->size); // gettimeofday(&fft_setup_time, NULL); // int fftw_thread = fftw_init_threads(); diff --git a/src/c_lib/sandbox/sandbox.pxd b/src/c_lib/sandbox/sandbox.pxd index 09dd85ee1..cd67c92fb 100644 --- a/src/c_lib/sandbox/sandbox.pxd +++ b/src/c_lib/sandbox/sandbox.pxd @@ -12,22 +12,22 @@ from libc.stdint cimport int32_t cdef extern from "schemes.h": ctypedef struct MRS_averaging_scheme: - unsigned int total_orientations - unsigned int integration_density - unsigned int integration_volume + int total_orientations + int integration_density + int integration_volume MRS_averaging_scheme * MRS_create_averaging_scheme( - unsigned int integration_density, + int integration_density, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, + int n_gamma, + int integration_volume, bool_t interpolation) MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( double *alpha, double *beta, - double *weight, unsigned int n_angles, + double *weight, int n_angles, bool_t allow_4th_rank, - const unsigned int position_size, + const int position_size, int32_t *positions, bool_t interpolation) @@ -36,11 +36,11 @@ cdef extern from "schemes.h": cdef extern from "mrsimulator.h": ctypedef struct MRS_plan: - unsigned int number_of_sidebands + int number_of_sidebands double rotor_frequency_in_Hz double rotor_angle_in_rad - MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, + MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, bool_t allow_4th_rank) diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index e9281938f..a79241cab 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -65,11 +65,11 @@ cdef extern from "angular_momentum/wigner_matrix.h": void wigner_dm0_vector(const int l, const double beta, double *R_out) - void get_exp_Im_angle(const unsigned int octant_orientations, + void get_exp_Im_angle(const int octant_orientations, const bool_t allow_4th_rank, void *exp_Im_angle, double delta_alpha) - void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, double *wigner_2j_matrices, void *R2, + void __batch_wigner_rotation(const int octant_orientations, + const int n_octants, double *wigner_2j_matrices, void *R2, double *wigner_4j_matrices, void *R4, void *exp_Im_alpha, void *w2, void *w4) @@ -90,7 +90,7 @@ cdef extern from "interpolation.h": double amp, double *spec, int points, - unsigned int iso_intrp) + int iso_intrp) void triangle_interpolation1D_linear( double freq1, @@ -119,26 +119,26 @@ cdef extern from "interpolation.h": double *spec, int m0, int m1, - unsigned int iso_intrp) + int iso_intrp) void octahedronInterpolation( double *spec, double *freq, int nt, double *amp, - const unsigned int stride, + const int stride, int m) cdef extern from "mrsimulator.h": void get_sideband_phase_components( - unsigned int number_of_sidebands, + int number_of_sidebands, double spin_frequency, double *pre_phase) # ctypedef struct MRS_plan # MRS_plan *MRS_create_plan( -# unsigned int integration_density, +# int integration_density, # int number_of_sidebands, # double rotor_frequency_in_Hz, # double rotor_angle_in_rad, @@ -172,7 +172,7 @@ cdef extern from "method.h": double increment # Increment of coordinates along the dimension. double coordinates_offset # Start coordinate of the dimension. MRS_event *events # Holds a list of events. - unsigned int n_events # The number of events. + int n_events # The number of events. # MRS_dimension *MRS_create_dimensions( # MRS_averaging_scheme *scheme, @@ -182,7 +182,7 @@ cdef extern from "method.h": # double *magnetic_flux_density_in_T, # double *rotor_frequency_in_Hz, # double *rotor_angle_in_rad, - # unsigned int n_events, + # int n_events, # int number_of_sidebands) cdef extern from "simulation.h": @@ -201,14 +201,14 @@ cdef extern from "simulation.h": # second order quad Hamiltonian. # spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, + int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, # The transition as transition[0] = mi and transition[1] = mf double *transition, int integration_density, - unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere. + int integration_volume, # 0-octant, 1-hemisphere, 2-sphere. ) cdef extern from "frequency/spatial_orientation_tensor_components.h": diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index 40bce4126..2d7a3f690 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -111,7 +111,7 @@ def __wigner_rotation_2(int l, np.ndarray[double] cos_alpha, @cython.boundscheck(False) @cython.wraparound(False) def get_exp_Im_angle(int n, np.ndarray[double] cos_alpha, bool_t allow_4th_rank): - cdef unsigned int n_ = n + cdef int n_ = n cdef np.ndarray[double complex] exp_Im_angle = np.empty(4*n, dtype=np.complex128) exp_Im_angle[3*n:] = cos_alpha + 1j*np.sqrt(1.0 - cos_alpha**2) clib.get_exp_Im_angle(n_, allow_4th_rank, &exp_Im_angle[0], 0.0) @@ -120,7 +120,7 @@ def get_exp_Im_angle(int n, np.ndarray[double] cos_alpha, bool_t allow_4th_rank) @cython.boundscheck(False) @cython.wraparound(False) -def pre_phase_components(unsigned int number_of_sidebands, double rotor_frequency_in_Hz): +def pre_phase_components(int number_of_sidebands, double rotor_frequency_in_Hz): cdef int n1 = 4 * number_of_sidebands cdef np.ndarray[double] pre_phase = np.zeros(2*n1, dtype=np.float64) clib.get_sideband_phase_components(number_of_sidebands, rotor_frequency_in_Hz, &pre_phase[0]) @@ -157,7 +157,7 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): :return amp: The amplitude at the given $\alpha$ and $\beta$. """ nt = integration_density - cdef unsigned int octant_orientations = int((nt+1) * (nt+2)/2) + cdef int octant_orientations = int((nt+1) * (nt+2)/2) cdef bool_t interploation = True cdef np.ndarray[double complex] exp_I_alpha = np.empty(octant_orientations, dtype=np.complex128) @@ -171,7 +171,7 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): @cython.boundscheck(False) @cython.wraparound(False) -def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, const unsigned int stride=1): +def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, const int stride=1): cdef int i cdef int number_of_sidebands = int(amp.shape[0]) for i in range(number_of_sidebands): @@ -197,7 +197,7 @@ def triangle_interpolation1D(vector, np.ndarray[double, ndim=1] spectrum_amp, default value is 0. :ivar type: Linear or Gaussian interpolation for delta functions. """ - cdef int points = np.asarray([spectrum_amp.size/2], dtype=np.int32)[0] + cdef int points = np.asarray([spectrum_amp.size/2], dtype=np.intc)[0] cdef np.ndarray[double, ndim=1] f_vector = np.asarray(vector, dtype=np.float64) cdef double f1 = f_vector[0] @@ -222,7 +222,7 @@ def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectr :ivar vector2: 1-D array of three points. :ivar spectrum_amp: A numpy array of amplitudes. This array is the output. """ - shape = np.asarray([spectrum_amp.shape[0], spectrum_amp.shape[1]/2], dtype=np.int32) + shape = np.asarray([spectrum_amp.shape[0], spectrum_amp.shape[1]/2], dtype=np.intc) # cdef np.ndarray[int, ndim=1] points = shape cdef np.ndarray[double, ndim=1] f1_vector = np.asarray(vector1, dtype=np.float64) cdef np.ndarray[double, ndim=1] f2_vector = np.asarray(vector2, dtype=np.float64) @@ -241,8 +241,8 @@ def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectr @cython.boundscheck(False) @cython.wraparound(False) -def __batch_wigner_rotation(unsigned int octant_orientations, - unsigned int n_octants, +def __batch_wigner_rotation(int octant_orientations, + int n_octants, np.ndarray[double] wigner_2j_matrices, np.ndarray[double complex] R2, np.ndarray[double] wigner_4j_matrices, From 2b2259839963d7536f49c5cbcddc5a948c7ed6b8 Mon Sep 17 00:00:00 2001 From: deepanshs Date: Wed, 28 Aug 2024 19:28:58 -0400 Subject: [PATCH 24/27] change all uint to int --- src/c_lib/lib/frequency_averaging.c | 17 +++-- src/c_lib/lib/interpolation.c | 108 +++++++++++++--------------- src/c_lib/lib/mrsimulator.c | 59 ++++++++------- src/c_lib/lib/simulation.c | 20 +++--- 4 files changed, 97 insertions(+), 107 deletions(-) diff --git a/src/c_lib/lib/frequency_averaging.c b/src/c_lib/lib/frequency_averaging.c index 09a307e38..8c9b41906 100644 --- a/src/c_lib/lib/frequency_averaging.c +++ b/src/c_lib/lib/frequency_averaging.c @@ -64,10 +64,9 @@ static inline void sideband_amplitude(int npts, int n_octant, complex128 *a11, } void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, unsigned int iso_intrp, - complex128 *exp_I_phase) { - unsigned int i, j, k1, address, ptr, gamma_idx; - unsigned int nt = scheme->integration_density, npts = scheme->octant_orientations; + double *spec, int iso_intrp, complex128 *exp_I_phase) { + int i, j, k1, address, ptr, gamma_idx; + int nt = scheme->integration_density, npts = scheme->octant_orientations; double offset_0, offset; double *freq, *phase_ptr, *amps = dimensions->freq_amplitude; @@ -158,10 +157,10 @@ void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * } void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, double *affine_matrix, - unsigned int iso_intrp, complex128 *exp_I_phase) { - unsigned int i, k, j, address, gamma_idx; - unsigned int npts = scheme->octant_orientations, ptr; + double *spec, double *affine_matrix, int iso_intrp, + complex128 *exp_I_phase) { + int i, k, j, address, gamma_idx; + int npts = scheme->octant_orientations, ptr; int *fft1_index, *fft2_index, n_min, n_max; MRS_plan *planA, *planB; @@ -170,7 +169,7 @@ void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * double offset0, offset1, offsetA, offsetB; double *freq0, *freq1; double norm0, norm1; - complex128 *res_t = malloc_complex128(npts); + complex128 *res_t = malloc_complex128((size_t)npts); bool user_defined = scheme->user_defined, interpolation = scheme->interpolation; offset0 = dimensions[0].R0_offset; diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index e4c82600c..78f9f14de 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -23,7 +23,7 @@ * @param amp The area corresponding to the frequency coordinate. * @param spec1 A pointer to the intensity vector. */ -static void inline delta_fn_linear_interpolation(const double freq, const int points, +inline static void delta_fn_linear_interpolation(const double freq, const int points, double amp, double *spec1) { int p = (int)freq; if (p >= points || p < 0) return; @@ -52,7 +52,7 @@ static void inline delta_fn_linear_interpolation(const double freq, const int po } } -static void inline delta_fn_gauss_interpolation(const double freq, const int points, +inline static void delta_fn_gauss_interpolation(const double freq, const int points, double amp, double *spec) { double res, w, a0, a1, a2, a3, a4, sum, temp; int p = (int)(floor(freq)), index, ip2, ip1, im1, im2, pad = 2, p2 = 2 * p; @@ -261,7 +261,7 @@ static inline void __triangle_interpolation(double freq1, double freq2, double f } void triangle_interpolation1D(double freq1, double freq2, double freq3, double amp, - double *spec, int points, unsigned int iso_intrp) { + double *spec, int points, int iso_intrp) { if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { if (iso_intrp == 0) delta_fn_linear_interpolation(freq1, points, amp, spec); if (iso_intrp == 1) delta_fn_gauss_interpolation(freq1, points, amp, spec); @@ -357,8 +357,7 @@ void triangle_interpolation1D_gaussian(double freq1, double freq2, double freq3, // f10 ------------------ f11 bottom line static inline void quadrilateral_bin(double f00, double f11, double f10, double f01, double top, double bottom, double total, - double amp, double *spec, int m1, - unsigned int iso_intrp) { + double amp, double *spec, int m1, int iso_intrp) { double amp_bottom, amp_top, norm; if (total != 0) { @@ -376,8 +375,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, bool r_clip, bool r_clip_r, double top, double *f1, double *f2, double *x10, double *x11, int m1, - double *spec, - unsigned int iso_intrp) { + double *spec, int iso_intrp) { double f10, slope_diff, abs_sdiff, abs_sdiff_2, temp, df1, diff; double amp, denom, line_up, line_down; double x00, x01, f01_slope, f02_slope; @@ -479,8 +477,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, bool r_clip, bool r_clip_l, double top, double *f1, double *f2, double *x10, double *x11, int m1, - double *spec, - unsigned int iso_intrp) { + double *spec, int iso_intrp) { double f21, slope_diff, abs_sdiff, abs_sdiff_2, temp, df2, diff; double amp, denom, line_up, line_down; double x00, x01, f12_slope, f02_slope; @@ -579,12 +576,12 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, void triangle_interpolation2D(double freq11, double freq12, double freq13, double freq21, double freq22, double freq23, double amp, - double *spec, int m0, int m1, unsigned int iso_intrp) { + double *spec, int m0, int m1, int iso_intrp) { double top, t1, t2, diff, temp, n_i; int p, pmid, pmax, i, j; double freq10_01 = 0.0, freq11_02 = 0.0; - p = (int)(freq11); + p = (int)freq11; if (absd(freq11 - freq12) < TOL && absd(freq11 - freq13) < TOL) { if (p >= m0 || p < 0) return; @@ -813,11 +810,11 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // self.x0 + t1*xdelta, self.y0 + t1*ydelta) // return [clipped_line] -void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *amp, - const unsigned int stride, int n_spec, double *spec, - unsigned int iso_intrp) { - unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; +void octahedronDeltaInterpolation(const int nt, double freq, double *amp, + const int stride, int n_spec, double *spec, + int iso_intrp) { + int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address; local_index = nt - 1; @@ -844,10 +841,10 @@ void octahedronDeltaInterpolation(const unsigned int nt, double freq, double *am if (iso_intrp == 1) delta_fn_gauss_interpolation(freq, n_spec, amp1, spec); } -void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, const unsigned int stride, int m) { - unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; +void octahedronInterpolation(double *spec, double *freq, const int nt, double *amp, + const int stride, int m) { + int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq_address; /* Interpolate between 1d points by setting up triangles of unit area */ @@ -880,9 +877,8 @@ void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, } void generic_1d_triangle_interpolation(double *spec, double *freq, double *amp, int m, - const unsigned int position_size, - int32_t *positions) { - unsigned int pos_size = position_size; + const int position_size, int32_t *positions) { + int pos_size = position_size; int32_t p1, p2, p3; double amp_sum; @@ -896,24 +892,23 @@ void generic_1d_triangle_interpolation(double *spec, double *freq, double *amp, } } -void hist1d(double *spec, const unsigned int freq_size, double *freq, double *amp, - int m) { - unsigned int i = 0, ix; +void hist1d(double *spec, const int freq_size, double *freq, double *amp, int m) { + int i = 0, ix; double temp_freq; for (i = 0; i < freq_size; i++) { temp_freq = freq[i]; if (temp_freq >= 0 && temp_freq < m) { - ix = (unsigned int)temp_freq; + ix = (int)temp_freq; spec[2 * ix] += amp[i]; } } } -void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, - double *amp_real, double *amp_imag, int dimension_count, - const unsigned int position_size, int32_t *positions, - const unsigned int nt, bool user_defined, bool interpolation) { +void one_d_averaging(double *spec, const int freq_size, double *freq, double *amp_real, + double *amp_imag, int dimension_count, const int position_size, + int32_t *positions, const int nt, bool user_defined, + bool interpolation) { if (!user_defined) { if (interpolation) { octahedronInterpolation(spec, freq, nt, amp_real, 1, dimension_count); @@ -930,9 +925,9 @@ void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, } } -void generic_1d_triangle_average(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, - const unsigned int position_size, int32_t *positions) { +void generic_1d_triangle_average(double *spec, const int freq_size, double *freq, + double *amp, int m, const int position_size, + int32_t *positions) { if (positions == NULL) { hist1d(spec, freq_size, freq, amp, m); } else { @@ -940,11 +935,11 @@ void generic_1d_triangle_average(double *spec, const unsigned int freq_size, } } -void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, - double *freq2, double *amp, int amp_stride, - const unsigned int position_size, int32_t *positions, - int dimension0_count, int dimension1_count, unsigned int iso_intrp, - const unsigned int nt, bool user_defined, bool interpolation) { +void two_d_averaging(double *spec, const int freq_size, double *freq1, double *freq2, + double *amp, int amp_stride, const int position_size, + int32_t *positions, int dimension0_count, int dimension1_count, + int iso_intrp, const int nt, bool user_defined, + bool interpolation) { if (!user_defined) { if (interpolation) { // Perform tenting on every sideband order over all orientations @@ -963,10 +958,9 @@ void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, void generic_2d_triangle_interpolation(double *spec, double *freq1, double *freq2, double *amp, int amp_stride, - const unsigned int position_size, - int32_t *positions, int m0, int m1, - unsigned int iso_intrp) { - unsigned int pos_size = position_size; + const int position_size, int32_t *positions, + int m0, int m1, int iso_intrp) { + int pos_size = position_size; int32_t p1, p2, p3; double amp_sum; @@ -981,28 +975,27 @@ void generic_2d_triangle_interpolation(double *spec, double *freq1, double *freq } } -void hist2d(double *spec, const unsigned int freq_size, double *freq1, double *freq2, +void hist2d(double *spec, const int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, int m0, int m1) { - unsigned int i = 0, ix, iy, hist_index; + int i = 0, ix, iy, hist_index; double temp_freq_1, temp_freq_2; for (i = 0; i < freq_size; i++) { temp_freq_1 = freq1[i]; temp_freq_2 = freq2[i]; if (temp_freq_1 >= 0 && temp_freq_1 < m0 && temp_freq_2 >= 0 && temp_freq_2 < m1) { - ix = (unsigned int)temp_freq_1; - iy = (unsigned int)temp_freq_2; + ix = (int)temp_freq_1; + iy = (int)temp_freq_2; hist_index = iy + m1 * ix; spec[2 * hist_index] += amp[i * amp_stride]; } } } -void generic_2d_triangle_average(double *spec, const unsigned int freq_size, - double *freq1, double *freq2, double *amp, - int amp_stride, int m0, int m1, - const unsigned int position_size, int32_t *positions, - unsigned int iso_intrp) { +void generic_2d_triangle_average(double *spec, const int freq_size, double *freq1, + double *freq2, double *amp, int amp_stride, int m0, + int m1, const int position_size, int32_t *positions, + int iso_intrp) { if (positions == NULL) { hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1); } else { @@ -1011,12 +1004,11 @@ void generic_2d_triangle_average(double *spec, const unsigned int freq_size, } } -void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - const unsigned int nt, double *amp, - const unsigned int stride, int m0, int m1, - unsigned int iso_intrp) { - unsigned int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; +void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, const int nt, + double *amp, const int stride, int m0, int m1, + int iso_intrp) { + int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq1_address, *freq2_address; /* Interpolate between 1d points by setting up triangles of unit area */ diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index 707e4ed35..7d88aa7c2 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -81,8 +81,7 @@ void MRS_plan_release_temp_storage(MRS_plan *the_plan) { * 4) creating the fftw plan, 4) allocating buffer for storing the evaluated frequencies * and their respective amplitudes. */ -MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, - unsigned int number_of_sidebands, +MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, bool allow_4th_rank) { MRS_plan *plan = malloc(sizeof(MRS_plan)); @@ -113,7 +112,7 @@ MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, * Normalizing amplitudes from the spherical averaging scheme by the number of * sidebands square times the number of octants. */ - plan->norm_amplitudes = malloc_double(scheme->octant_orientations); + plan->norm_amplitudes = malloc_double((size_t)scheme->octant_orientations); cblas_dcopy(scheme->octant_orientations, scheme->amplitudes, 1, plan->norm_amplitudes, 1); double scale = (1.0 / (double)(plan->number_of_sidebands * plan->number_of_sidebands * @@ -135,7 +134,7 @@ MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, */ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, double rotor_frequency_in_Hz) { - unsigned int size_4; + int size_4; plan->rotor_frequency_in_Hz = rotor_frequency_in_Hz; plan->is_static = (rotor_frequency_in_Hz < 1.0e-3) ? true : false; @@ -150,7 +149,7 @@ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, * @see get_sideband_phase_components() */ size_4 = 4 * plan->number_of_sidebands; - plan->pre_phase = malloc_complex128(size_4); + plan->pre_phase = malloc_complex128((size_t)size_4); get_sideband_phase_components(plan->number_of_sidebands, rotor_frequency_in_Hz, (double *)plan->pre_phase); } @@ -164,14 +163,14 @@ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, */ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_in_rad, bool allow_4th_rank) { - unsigned int size_2, size_4, i, j; + int size_2, size_4, i, j; plan->rotor_angle_in_rad = rotor_angle_in_rad; /** * Calculate wigner-2j d^2_{m,0} vector where m ∈ [-2, 2]. This vector is used to * rotate the second-rank tensors from the rotor frame to the lab frame. * @see wigner_dm0_vector() */ - plan->wigner_d2m0_vector = malloc_double(5); + plan->wigner_d2m0_vector = malloc_double((size_t)5); wigner_dm0_vector(2, rotor_angle_in_rad, plan->wigner_d2m0_vector); plan->wigner_d4m0_vector = NULL; @@ -181,13 +180,13 @@ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_ * rotate the fourth-rank tensors from the rotor frame to the lab frame. * @see wigner_dm0_vector() */ - plan->wigner_d4m0_vector = malloc_double(9); + plan->wigner_d4m0_vector = malloc_double((size_t)9); wigner_dm0_vector(4, rotor_angle_in_rad, plan->wigner_d4m0_vector); } // pre_phase_2 is only calculated for m=-2 and -1 for l=2 rank tensor calculation. size_2 = 2 * plan->number_of_sidebands; - plan->pre_phase_2 = malloc_complex128(size_2); + plan->pre_phase_2 = malloc_complex128((size_t)size_2); /* Copy the pre_phase[m=-2 to 2] to pre_phase2 */ cblas_zcopy(size_2, (double *)(plan->pre_phase[2 * plan->number_of_sidebands]), 1, @@ -215,7 +214,7 @@ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_ /* pre_phase_4 is only calculated for m=-4, -3, -2, and -1 for l=4 rank tensor * calculation. */ size_4 = 4 * plan->number_of_sidebands; - plan->pre_phase_4 = malloc_complex128(size_4); + plan->pre_phase_4 = malloc_complex128((size_t)size_4); /* Copy the pre_phase[m=-4 to 4] to pre_phase4 */ cblas_zcopy(size_4, (double *)(plan->pre_phase), 1, (double *)(plan->pre_phase_4), 1); @@ -405,8 +404,8 @@ void MRS_get_normalized_frequencies_from_plan(MRS_averaging_scheme *scheme, double fraction, unsigned char is_spectral, double duration) { - unsigned int i, g_idx, ptr; - unsigned int freq_size = scheme->n_gamma * scheme->total_orientations; + int i, g_idx, ptr; + int freq_size = scheme->n_gamma * scheme->total_orientations; double temp, fraction_duration; double *f_complex, *local_frequency; @@ -441,7 +440,7 @@ void MRS_get_normalized_frequencies_from_plan(MRS_averaging_scheme *scheme, /* Normalized local anisotropic frequency contributions from the 2nd-rank tensor. */ for (g_idx = 0; g_idx < scheme->n_gamma; g_idx++) { ptr = scheme->total_orientations * g_idx; - temp = 2 * fraction_duration; + temp = 2.0 * fraction_duration; if (plan->is_static) { for (i = 0; i < 2; i++) { @@ -507,7 +506,7 @@ static inline void MRS_rotate_single_site_interaction_components( double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib // The pointer to freq contribs boolean. ) { - unsigned int i, n_sites = sites->number_of_sites; + int i, n_sites = sites->number_of_sites; double larmor_freq_in_MHz, larmor_freq_in_Hz; float *mf = &transition[n_sites], *mi = transition; double F2_shield[10]; @@ -586,14 +585,14 @@ static inline void quad_coupling_cross_terms( coupling_struct *couplings, // Pointer to a list of couplings within a spin system. site_struct *sites, // Pointer to a list of sites within a spin system. int site_index, // site index - unsigned int coupling_index, // coupled index - double *F0, // The F0 components. - complex128 *F2, // The F2 components. - complex128 *F4, // The F4 components. - double *F0_temp, // The temporary F0 components. - complex128 *F2_temp, // The temporary F2 components. - complex128 *F4_temp, // The temporary F4 components. - double B0_in_T, // Magnetic flux density in T. + int coupling_index, // coupled index + double *F0, // The F0 components. + complex128 *F2, // The F2 components. + complex128 *F4, // The F4 components. + double *F0_temp, // The temporary F0 components. + complex128 *F2_temp, // The temporary F2 components. + complex128 *F4_temp, // The temporary F4 components. + double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib, // The pointer to freq contribs boolean. float mIi, // inital transition state of site coupled to the quad site (p) float mSqi, // inital transition state of the quad site (d) @@ -640,7 +639,7 @@ static inline void MRS_rotate_coupled_site_interaction_components( coupling_struct *couplings, // Pointer to a list of couplings within a spin system. site_struct *sites, // Pointer to a list of sites within a spin system. float *transition, // The spin transition. - unsigned int n_sites, // The number of sites. + int n_sites, // The number of sites. double *F0, // The F0 components. complex128 *F2, // The F2 components. complex128 *F4, // The F4 components. @@ -650,7 +649,7 @@ static inline void MRS_rotate_coupled_site_interaction_components( double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib // The pointer to freq contribs boolean. ) { - unsigned int i, j = 0, n_couplings = couplings->number_of_couplings; + int i, j = 0, n_couplings = couplings->number_of_couplings; int site_index_I, site_index_S; float mIf, mSf, mIi, mSi; @@ -757,15 +756,15 @@ void MRS_rotate_components_from_PAS_to_common_frame( * = scale [-[cos(m ωr t) -1] +Isin(m ωr t)] * That is, pre_phase[-m] = -Re(pre_phase[m]) + Im(pre_phase[m]) */ -void get_sideband_phase_components_2(unsigned int number_of_sidebands, +void get_sideband_phase_components_2(int number_of_sidebands, double rotor_frequency_in_Hz, complex128 *pre_phase) { int m, i; double spin_angular_freq, tau, scale; - double *input = malloc_double(number_of_sidebands); - double *ones = malloc_double(number_of_sidebands); - double *phase = malloc_double(number_of_sidebands); + double *input = malloc_double((size_t)number_of_sidebands); + double *ones = malloc_double((size_t)number_of_sidebands); + double *phase = malloc_double((size_t)number_of_sidebands); vm_double_ones(number_of_sidebands, ones); vm_double_arange(number_of_sidebands, input); @@ -833,11 +832,11 @@ void get_sideband_phase_components_2(unsigned int number_of_sidebands, * as the leading dimension. The first number_of_sidebands entries corresponds to * m_wr=-4. */ -void get_sideband_phase_components(unsigned int number_of_sidebands, +void get_sideband_phase_components(int number_of_sidebands, double sample_rotation_frequency, double *restrict pre_phase) { double spin_angular_freq, tau, wrt, pht, scale; - unsigned int step, m; + int step, m; // Calculate the spin angular frequency spin_angular_freq = sample_rotation_frequency * CONST_2PI; diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index 323926214..fb9e46ad1 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -50,7 +50,7 @@ void __mrsimulator_core( MRS_dimension *dimensions, // Pointer to MRS_dimension structure. MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - unsigned int iso_intrp, // Isotropic interpolation scheme (linear | Gaussian) + int iso_intrp, // Isotropic interpolation scheme (linear | Gaussian) unsigned char *freq_contrib, // A list of freq_contrib booleans. double *affine_matrix // Affine transformation matrix. ) { @@ -70,22 +70,22 @@ void __mrsimulator_core( // sizeof(transition_pathway_weight[0])); printf("n_dimension (int) %zu\n\n", // sizeof(n_dimension)); - // printf("iso_intrp (unsigned int) %zu\n", sizeof(iso_intrp)); + // printf("iso_intrp (int) %zu\n", sizeof(iso_intrp)); // printf("freq_contrib (unsigned char) %zu\n", sizeof(freq_contrib[0])); // printf("affine_matrix double %zu\n\n", sizeof(affine_matrix[0])); - // printf("n_sites site (unsigned int) %zu\n", sizeof(sites->number_of_sites)); + // printf("n_sites site (int) %zu\n", sizeof(sites->number_of_sites)); // printf("n_sites spin (float) %zu\n", sizeof(sites->spin[0])); // printf("n_sites gyromagnetic_ratio (double) %zu\n\n", // sizeof(sites->gyromagnetic_ratio[0])); - // printf("n_coupling n_c (unsigned int) %zu\n", + // printf("n_coupling n_c (int) %zu\n", // sizeof(couplings->number_of_couplings)); printf("n_coupling site_index (int) // %zu\n", sizeof(couplings->site_index[0])); printf("isotropic_j_in_Hz (double) // %zu\n\n", sizeof(couplings->isotropic_j_in_Hz[0])); unsigned char is_spectral; - unsigned int evt; + int evt; int dim, total_pts = scheme->n_gamma * scheme->total_orientations; double B0_in_T, fraction, duration; vm_double_zeros(total_pts, scheme->phase); @@ -192,16 +192,16 @@ void mrsimulator_core( int quad_second_order, // Quad theory for second order, // spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, // The number of sidebands - double rotor_frequency_in_Hz, // The rotor spin frequency - double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis + int number_of_sidebands, // The number of sidebands + double rotor_frequency_in_Hz, // The rotor spin frequency + double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis // Pointer to the a list of transitions. float *transition_pathway, double *transition_pathway_weight, // powder orientation average int integration_density, // The number of triangle along the edge of octahedron - unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. - unsigned int interpolate_type, unsigned char *freq_contrib, double *affine_matrix) { + int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix) { // int num_process = openblas_get_num_procs(); // int num_threads = openblas_get_num_threads(); // openblas_set_num_threads(1); From a47715eac2b1e07eb277c32fc1bfdca778bfd71b Mon Sep 17 00:00:00 2001 From: deepanshs Date: Wed, 28 Aug 2024 20:30:29 -0400 Subject: [PATCH 25/27] fix --- src/c_lib/lib/simulation.c | 26 +------------------------- 1 file changed, 1 insertion(+), 25 deletions(-) diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index fb9e46ad1..a36096087 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -59,30 +59,6 @@ void __mrsimulator_core( et al. `Computation of Orientational Averages in Solid-State NMR by Gaussian Spherical Quadrature` JMR, 132, 1998. https://doi.org/10.1006/jmre.1998.1427 */ - // printf("int %zu\n", sizeof(int)); - // printf("double %zu\n", sizeof(double)); - // printf("float %zu\n", sizeof(float)); - // printf("unsigned int %zu\n", sizeof(unsigned int)); - // printf("unsigned char %zu\n\n", sizeof(unsigned char)); - - // printf("transition_pathway (float) %zu\n", sizeof(transition_pathway[0])); - // printf("transition_pathway_weight (double) %zu\n", - // sizeof(transition_pathway_weight[0])); printf("n_dimension (int) %zu\n\n", - // sizeof(n_dimension)); - - // printf("iso_intrp (int) %zu\n", sizeof(iso_intrp)); - // printf("freq_contrib (unsigned char) %zu\n", sizeof(freq_contrib[0])); - // printf("affine_matrix double %zu\n\n", sizeof(affine_matrix[0])); - - // printf("n_sites site (int) %zu\n", sizeof(sites->number_of_sites)); - // printf("n_sites spin (float) %zu\n", sizeof(sites->spin[0])); - // printf("n_sites gyromagnetic_ratio (double) %zu\n\n", - // sizeof(sites->gyromagnetic_ratio[0])); - - // printf("n_coupling n_c (int) %zu\n", - // sizeof(couplings->number_of_couplings)); printf("n_coupling site_index (int) - // %zu\n", sizeof(couplings->site_index[0])); printf("isotropic_j_in_Hz (double) - // %zu\n\n", sizeof(couplings->isotropic_j_in_Hz[0])); unsigned char is_spectral; int evt; @@ -107,7 +83,7 @@ void __mrsimulator_core( // Loop over the dimension. for (dim = 0; dim < n_dimension; dim++) { // Reset the freqs to zero at the start of each spectral dimension. - cblas_dscal(total_pts, 0.0, dimensions[dim].local_frequency, 1); + vm_double_zeros(total_pts, dimensions[dim].local_frequency); dimensions[dim].R0_offset = 0.0; plan = dimensions[dim].events->plan; From 1f939d49137d97b445d1c64faecee794f3be2302 Mon Sep 17 00:00:00 2001 From: deepanshs <21365911+deepanshs@users.noreply.github.com> Date: Wed, 28 Aug 2024 20:40:19 -0400 Subject: [PATCH 26/27] cleanup --- .../workflows/continuous-integration-pip.yml | 145 +++++++++--------- setup.py | 1 + src/c_lib/base/base_model.pyx | 29 ---- src/mrsimulator/simulator/__init__.py | 7 - 4 files changed, 73 insertions(+), 109 deletions(-) diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index 9b4e324cc..b8c50c0c5 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -38,49 +38,48 @@ jobs: # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide flake8 . --count --exit-zero --max-complexity=8 --max-line-length=88 --statistics --exclude="examples/* *.npy docs/* *.pyx *.pxd" - # testing_unix: - # needs: [code_lint] - # runs-on: "ubuntu-latest" - # strategy: - # matrix: - # python-version: [3.9, "3.10", "3.11", "3.12"] - - # steps: - # - name: Checkout - # uses: actions/checkout@v4 - - # - name: Set up Python ${{ matrix.python-version }} - # uses: actions/setup-python@v5 - # with: - # python-version: ${{ matrix.python-version }} - - # - name: Install linux system dependencies - # run: | - # sudo apt-get install --yes libopenblas-dev libfftw3-dev - # gcc -v - - # - 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@v4.5.0 + testing_unix: + needs: [code_lint] + runs-on: "ubuntu-latest" + strategy: + matrix: + python-version: [3.9, "3.10", "3.11", "3.12"] + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Install linux system dependencies + run: | + sudo apt-get install --yes libopenblas-dev libfftw3-dev + gcc -v + + - 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@v4.5.0 testing_mac_arm: needs: [code_lint] runs-on: "macos-latest" strategy: matrix: - python-version: ["3.10"] - # python-version: [3.9, "3.10", "3.11", "3.12"] + python-version: [3.9, "3.10", "3.11", "3.12"] steps: - name: Checkout @@ -141,7 +140,7 @@ jobs: pip install -r requirements-dev.txt - name: Build and install package from source - run: CC=gcc-12 python setup.py develop + run: python setup.py develop - name: Test with pytest run: pytest --cov=./ --cov-report=xml @@ -149,38 +148,38 @@ jobs: - name: Upload coverage uses: codecov/codecov-action@v4.5.0 - # testing_windows: - # needs: [code_lint] - # runs-on: "windows-latest" - # strategy: - # matrix: - # python-version: [3.9, "3.10", "3.11", "3.12"] - - # steps: - # - name: Checkout - # uses: actions/checkout@v4 - # - name: Setup Miniconda - # uses: conda-incubator/setup-miniconda@v3 - # env: - # ACTIONS_ALLOW_UNSECURE_COMMANDS: "true" - # with: - # auto-update-conda: true - # auto-activate-base: false - # miniconda-version: "latest" - # python-version: ${{ matrix.python-version }} - # environment-file: environment-dev.yml - # activate-environment: mrsimulator-dev - # - run: | - # conda --version - # which python - # - name: Build and install package from source - # shell: pwsh - # run: | - # conda --version - # which python - # python setup.py develop - # - name: Test with pytest - # shell: pwsh - # run: pytest --cov=./ --cov-report=xml - # - name: Upload coverage - # uses: codecov/codecov-action@v4.5.0 + testing_windows: + needs: [code_lint] + runs-on: "windows-latest" + strategy: + matrix: + python-version: [3.9, "3.10", "3.11", "3.12"] + + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Setup Miniconda + uses: conda-incubator/setup-miniconda@v3 + env: + ACTIONS_ALLOW_UNSECURE_COMMANDS: "true" + with: + auto-update-conda: true + auto-activate-base: false + miniconda-version: "latest" + python-version: ${{ matrix.python-version }} + environment-file: environment-dev.yml + activate-environment: mrsimulator-dev + - run: | + conda --version + which python + - name: Build and install package from source + shell: pwsh + run: | + conda --version + which python + python setup.py develop + - name: Test with pytest + shell: pwsh + run: pytest --cov=./ --cov-report=xml + - name: Upload coverage + uses: codecov/codecov-action@v4.5.0 diff --git a/setup.py b/setup.py index a1734898f..0a5cdef8c 100644 --- a/setup.py +++ b/setup.py @@ -192,6 +192,7 @@ def __init__(self): "-fcommon", "-Wall", "-Wextra", + "-Wconversion", ] self.extra_link_args += ["-lm"] diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index 7fed0ff87..3833717b4 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -447,35 +447,6 @@ def core_simulator(method, # print('pathway', transition_pathway_c) # print('weight', transition_pathway_weight_c) # print('pathway_count, inc', pathway_count, pathway_increment) - - if debug: - print("int ", sizeof(int)); - print("double ", sizeof(double)); - print("float ", sizeof(float)); - print("unsigned int ", sizeof(unsigned int)); - print("unsigned char ", sizeof(unsigned char)); - print() - - print("transition_pathway (float) ", sizeof(transition_pathway_c[0])); - print("transition_pathway_weight (double) ", sizeof(transition_pathway_weight_c[0])); - print("n_dimension (int) ", sizeof(n_dimension)); - print() - - print("iso_intrp (unsigned int) ", sizeof(isotropic_interpolation)); - print("freq_contrib (unsigned char) ", sizeof(f_contrib[0])); - print("affine_matrix (double) ", sizeof(affine_matrix_c[0])); - print() - - print("n_sites site (unsigned int) ", sizeof(sites_c.number_of_sites)); - print("n_sites spin (float) ", sizeof(sites_c.spin[0])); - print("n_sites gyromagnetic_ratio (double) ", sizeof(sites_c.gyromagnetic_ratio[0])); - print() - - print("n_coupling n_c (unsigned int) ", sizeof(couplings_c.number_of_couplings)); - print("n_coupling site_index (int) ", sizeof(couplings_c.site_index[0])); - print("isotropic_j_in_Hz (double) ", sizeof(couplings_c.isotropic_j_in_Hz[0])); - print() - for trans__ in range(pathway_count): clib.__mrsimulator_core( &[0], # as complex array diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index aad29de34..9cb2d8872 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -587,13 +587,6 @@ def _package_amp_after_simulation(self, method, amp, pack_as_csdm): bool pack_as_csdm: Packages the simulated spectrum as a CSDM object if true. Otherwise kept as a numpy array. """ - # if isinstance(amp[0], list): - # simulated_dataset = [] - # for item in amp: - # simulated_dataset += item - # if isinstance(amp[0], np.ndarray): - # simulated_dataset = [np.asarray(amp).sum(axis=0)] - method.simulation = ( self._as_csdm_object(amp, method) if pack_as_csdm else np.asarray(amp) ) From d07eb3ffe8b004d5d615af12907f7c96ff6502e3 Mon Sep 17 00:00:00 2001 From: unknown <21365911+deepanshs@users.noreply.github.com> Date: Thu, 29 Aug 2024 23:52:06 -0400 Subject: [PATCH 27/27] clean up --- setup.py | 8 ++++---- src/c_lib/lib/mrsimulator.c | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/setup.py b/setup.py index 0a5cdef8c..5ff8cebb6 100644 --- a/setup.py +++ b/setup.py @@ -66,7 +66,7 @@ def check_if_header_exists(self, header): def conda_setup_for_windows(self): self.libraries += ["fftw3", "openblas"] - self.extra_compile_args = ["/DUSE_OPENBLAS"] + self.extra_compile_args += ["/DUSE_OPENBLAS"] print(sys.version) loc = dirname(sys.executable) @@ -93,7 +93,7 @@ def conda_setup_for_unix(self): self.include_dirs += self.check_valid_path([join(loc, "include")]) self.library_dirs += self.check_valid_path([join(loc, "lib")]) - self.extra_compile_args = ["-O3", "-ffast-math", "-DUSE_OPENBLAS"] + self.extra_compile_args += ["-O3", "-ffast-math", "-DUSE_OPENBLAS"] self.libraries += ["fftw3", "openblas"] def on_exit_message(self, blas_lib, fftw_lib): @@ -132,8 +132,8 @@ class WindowsSetup(Setup): def __init__(self): super().__init__() - self.extra_link_args += ["-Wl"] - self.extra_compile_args = ["-DFFTW_DLL"] + self.extra_link_args += [] + self.extra_compile_args += ["-DFFTW_DLL", "/permissive-"] # if use_mkl: # self.mkl_blas_info() diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index 7d88aa7c2..c3388f280 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -414,7 +414,7 @@ void MRS_get_normalized_frequencies_from_plan(MRS_averaging_scheme *scheme, fraction_duration = dim->inverse_increment * fraction; } else { local_frequency = dim->local_phase; - cblas_dscal(freq_size, 0.0, local_frequency, 1); + vm_double_zeros(freq_size, local_frequency); fraction_duration = CONST_2PI * duration; }