diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index de05545..992d3ce 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -12,8 +12,8 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-13, macos-14, ubuntu-20.04, ubuntu-22.04, windows-2019, windows-2022] - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + os: [macos-14, macos-15, macos-26, ubuntu-22.04, ubuntu-24.04, windows-2022, windows-2025] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] include: - os: macos-latest env: CC="clang" CXX="clang++" @@ -23,22 +23,20 @@ jobs: env: CC="mingw64-gcc" CXX="mingw64-g++" steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v5 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v2 + uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} - name: Install dependencies run: | - python -m pip install --upgrade pip - python -m pip install --upgrade setuptools wheel + python -m pip install --upgrade pip setuptools wheel python -m pip install --upgrade flake8 pytest pytest-cov coverage spectres python -m pip install --upgrade -r requirements.txt - python setup.py build python -m pip install -e . - name: Test with pytest run: | coverage run -m pytest --disable-warnings -p no:warnings test - name: Coveralls - if: startsWith(matrix.os, 'ubuntu') && endsWith(matrix.python-version, '3.10') + if: startsWith(matrix.os, 'ubuntu') && endsWith(matrix.python-version, '3.14') uses: coverallsapp/github-action@v2 diff --git a/.github/workflows/cbuild.yml b/.github/workflows/cbuild.yml index 1208202..553faf3 100644 --- a/.github/workflows/cbuild.yml +++ b/.github/workflows/cbuild.yml @@ -11,31 +11,219 @@ on: jobs: build_wheels: - name: Build wheels on ${{ matrix.os }} + name: Build wheel for ${{ matrix.os }}-${{ matrix.cibw_python }}-${{ matrix.cibw_arch }} runs-on: ${{ matrix.os }} strategy: + fail-fast: false # Continue other builds even if one fails matrix: - os: [ubuntu-latest, windows-latest, macos-13, macos-14] + include: + # Linux x86_64 + - os: ubuntu-latest + cibw_python: cp39 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp310 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp311 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp312 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp313 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp314 + cibw_arch: x86_64 + # macOS ARM64 (Apple Silicon) + - os: macos-14 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp314 + cibw_arch: arm64 + + - os: macos-15 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp314 + cibw_arch: arm64 + + - os: macos-26 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp314 + cibw_arch: arm64 + + # macOS Intel x86_64 + - os: macos-14 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp314 + cibw_arch: x86_64 + + - os: macos-15 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp314 + cibw_arch: x86_64 + + - os: macos-26 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp314 + cibw_arch: x86_64 + + # Windows AMD64 + - os: windows-latest + cibw_python: cp39 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp310 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp311 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp312 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp313 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp314 + cibw_arch: AMD64 + + # Windows x86 + - os: windows-latest + cibw_python: cp39 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp310 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp311 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp312 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp313 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp314 + cibw_arch: x86 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.x' - name: Build wheels - uses: pypa/cibuildwheel@v2.19.1 + uses: pypa/cibuildwheel@v3.2.1 env: - # Disable building PyPy wheels on all platforms + CIBW_BUILD: ${{ matrix.cibw_python }}-* + CIBW_ARCHS: ${{ matrix.cibw_arch }} CIBW_SKIP: pp* + CIBW_BUILD_VERBOSITY: 3 + + - name: Verify wheel + run: | + python -m pip install --upgrade pip twine + python -m twine check wheelhouse/*.whl + shell: bash - uses: actions/upload-artifact@v4 with: - name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }} + name: cibw-wheel-${{ matrix.os }}-${{ matrix.cibw_python }}-${{ matrix.cibw_arch }} path: ./wheelhouse/*.whl build_sdist: name: Build source distribution runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - name: Build sdist run: pipx run build --sdist @@ -48,7 +236,6 @@ jobs: upload_pypi: needs: [build_wheels, build_sdist] runs-on: ubuntu-latest - environment: pypi permissions: id-token: write # or, alternatively, upload to PyPI on every tag starting with 'v' (remove on: release above to use this) @@ -56,14 +243,19 @@ jobs: steps: - uses: actions/download-artifact@v4 with: - # unpacks all CIBW artifacts into dist/ pattern: cibw-* - path: dist - merge-multiple: true + path: artifacts + merge-multiple: false + + - name: Move all wheels and tarballs into dist/ + run: | + mkdir -p dist + find artifacts -type f -name "*.whl" -exec mv {} dist/ \; + find artifacts -type f -name "*.tar.gz" -exec mv {} dist/ \; + ls -R dist - uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ - password: ${{ secrets.test_pypi_password }} - repository_url: https://test.pypi.org/legacy/ - skip-existing: true \ No newline at end of file + repository-url: https://test.pypi.org/legacy/ + skip-existing: true + packages-dir: dist \ No newline at end of file diff --git a/.github/workflows/cibuild-pypi.yml b/.github/workflows/cibuild-pypi.yml index b7e2043..68cc3b6 100644 --- a/.github/workflows/cibuild-pypi.yml +++ b/.github/workflows/cibuild-pypi.yml @@ -12,31 +12,219 @@ on: jobs: build_wheels: - name: Build wheels on ${{ matrix.os }} + name: Build wheel for ${{ matrix.os }}-${{ matrix.cibw_python }}-${{ matrix.cibw_arch }} runs-on: ${{ matrix.os }} strategy: + fail-fast: false # Continue other builds even if one fails matrix: - os: [ubuntu-latest, windows-latest, macos-13, macos-14] + include: + # Linux x86_64 + - os: ubuntu-latest + cibw_python: cp39 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp310 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp311 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp312 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp313 + cibw_arch: x86_64 + - os: ubuntu-latest + cibw_python: cp314 + cibw_arch: x86_64 + # macOS ARM64 (Apple Silicon) + - os: macos-14 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-14 + cibw_python: cp314 + cibw_arch: arm64 + + - os: macos-15 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-15 + cibw_python: cp314 + cibw_arch: arm64 + + - os: macos-26 + cibw_python: cp39 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp310 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp311 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp312 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp313 + cibw_arch: arm64 + - os: macos-26 + cibw_python: cp314 + cibw_arch: arm64 + + # macOS Intel x86_64 + - os: macos-14 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-14 + cibw_python: cp314 + cibw_arch: x86_64 + + - os: macos-15 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-15 + cibw_python: cp314 + cibw_arch: x86_64 + + - os: macos-26 + cibw_python: cp39 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp310 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp311 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp312 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp313 + cibw_arch: x86_64 + - os: macos-26 + cibw_python: cp314 + cibw_arch: x86_64 + + # Windows AMD64 + - os: windows-latest + cibw_python: cp39 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp310 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp311 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp312 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp313 + cibw_arch: AMD64 + - os: windows-latest + cibw_python: cp314 + cibw_arch: AMD64 + + # Windows x86 + - os: windows-latest + cibw_python: cp39 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp310 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp311 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp312 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp313 + cibw_arch: x86 + - os: windows-latest + cibw_python: cp314 + cibw_arch: x86 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.x' - name: Build wheels - uses: pypa/cibuildwheel@v2.19.1 + uses: pypa/cibuildwheel@v3.2.1 env: - # Disable building PyPy wheels on all platforms + CIBW_BUILD: ${{ matrix.cibw_python }}-* + CIBW_ARCHS: ${{ matrix.cibw_arch }} CIBW_SKIP: pp* + CIBW_BUILD_VERBOSITY: 3 + + - name: Verify wheel + run: | + python -m pip install --upgrade pip twine + python -m twine check wheelhouse/*.whl + shell: bash - uses: actions/upload-artifact@v4 with: - name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }} + name: cibw-wheel-${{ matrix.os }}-${{ matrix.cibw_python }}-${{ matrix.cibw_arch }} path: ./wheelhouse/*.whl build_sdist: name: Build source distribution runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - name: Build sdist run: pipx run build --sdist @@ -49,7 +237,6 @@ jobs: upload_pypi: needs: [build_wheels, build_sdist] runs-on: ubuntu-latest - environment: pypi permissions: id-token: write if: github.event_name == 'release' && github.event.action == 'published' @@ -58,13 +245,18 @@ jobs: steps: - uses: actions/download-artifact@v4 with: - # unpacks all CIBW artifacts into dist/ pattern: cibw-* - path: dist - merge-multiple: true + path: artifacts + merge-multiple: false + + - name: Move all wheels and tarballs into dist/ + run: | + mkdir -p dist + find artifacts -type f -name "*.whl" -exec mv {} dist/ \; + find artifacts -type f -name "*.tar.gz" -exec mv {} dist/ \; + ls -R dist - uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ - password: ${{ secrets.pypi_password }} - repository_url: https://upload.pypi.org/legacy/ \ No newline at end of file + repository_url: https://upload.pypi.org/legacy/ + package_dir: dist \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..cee16cd --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,10 @@ +CHANGELOG +========= + +Changes in v1.1.0 +----------------- + +- Enhanced C implementation to handle multi-dimensional flux arrays (2D, 3D) +- Added test coverage for multi-dimensional flux arrays +- Updated Python version support (removed 3.8, added 3.13-3.14) +- Updated GitHub Actions versions and macOS runner configurations diff --git a/README.md b/README.md index 5ed6871..460db1b 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ This is a Python package written with C extension to provide significant perform ![alt text](https://github.com/cylammarco/SpectResC/blob/main/speed_test/speed_test.png?raw=true) -We keep the implementation as close to SpectRes as possible. As of SpectRes v2.2.0, we do not see discrepant results between using SpectRes and SpectResC. +We keep the implementation as close to SpectRes as possible. The speed test was performed with a single core of Intel(R) Xeon(R) Gold 6226R CPU @ 2.90GHz (Oct 2025) against SpectRes v2.2.2. ## Installation diff --git a/pyproject.toml b/pyproject.toml index d1c16ab..11258ee 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -2,6 +2,36 @@ requires = ["setuptools>=42", "wheel", "cffi>=1.0", "numpy"] build-backend = "setuptools.build_meta" +[project] +name = "spectresc" +version = "1.1.0.dev0" +description = "SpectRes in C" +readme = "README.md" +requires-python = ">=3.9" +license = {text = "GPL-3.0"} +authors = [{name = "Marco C Lam", email = "mlam@roe.ac.uk"}] +maintainers = [{name = "Marco C Lam", email = "mlam@roe.ac.uk"}] +dependencies = [ + "numpy", +] +classifiers = [ + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Programming Language :: C", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", + "Operating System :: Microsoft :: Windows", + "Operating System :: Unix", + "Operating System :: MacOS", +] + +[project.urls] +Download = "https://github.com/cylammarco/SpectResC" + [toolset.cffi] extra_compile_args = ["-std=c99"] extra_link_args = [] @@ -28,3 +58,7 @@ ldflags = [] [tool.coverage.run] relative_files = true + +[tool.setuptools] +packages = ["spectresc"] +package-dir = {"" = "src"} diff --git a/setup.py b/setup.py index 8cdaac9..2f7ca1c 100644 --- a/setup.py +++ b/setup.py @@ -1,44 +1,27 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -from pathlib import Path + +from setuptools import Extension, setup, find_packages +import sys import numpy -from setuptools import Extension, setup + +# Platform-specific compile args +if sys.platform == "win32": + extra_compile_args = ["/O2"] # Optimize for Windows +else: + extra_compile_args = ["-O3", "-fPIC"] # Optimize for Unix/macOS setup( name="spectresc", - maintainer="Marco C Lam", - maintainer_email="mlam@roe.ac.uk", - version="1.0.4", - install_requires=[ - "numpy", - ], - description="SpectRes in C", - long_description=Path("README.md").read_text(encoding="utf-8"), - long_description_content_type="text/markdown", - author="Marco C Lam", - download_url="https://github.com/cylammarco/SpectResC", - classifiers=[ - "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", - "Programming Language :: C", - "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", - "Programming Language :: Python :: 3.10", - "Programming Language :: Python :: 3.11", - "Programming Language :: Python :: 3.12", - "Operating System :: Microsoft :: Windows", - "Operating System :: Unix", - "Operating System :: MacOS", - ], - platforms=["Windows", "Linux", "Mac OS-X"], - zip_safe=False, - include_dirs=["/usr/local/lib", numpy.get_include()], + include_dirs=[numpy.get_include()], + package_dir={"": "src"}, # Tell setuptools packages are under src/ + packages=find_packages(where="src"), # Find packages in src/ ext_modules=[ Extension( - "spectresc", + "spectresc.spectresc", sources=["src/spectresc/spectres.c"], - extra_compile_args=["-O3", "-fPIC", "-shared"], + extra_compile_args=extra_compile_args, ) ], ) diff --git a/speed_test/plot_speed_test_results.py b/speed_test/plot_speed_test_results.py index 5310915..87168b7 100644 --- a/speed_test/plot_speed_test_results.py +++ b/speed_test/plot_speed_test_results.py @@ -27,7 +27,6 @@ axs[0].plot( insizes, p / n, - label=f"numba with output size: {size}", color=red_gradient[j], ) @@ -35,7 +34,6 @@ axs[0].plot( insizes, p / c, - label=f"C extension with output size: {size}", color=blue_gradient[j], ) @@ -50,6 +48,7 @@ insizes, p_e / n_e, color=red_gradient[k], + label=f"numba with output size: {size}", ) for k, (p_e, c_e, size) in enumerate( @@ -63,12 +62,13 @@ insizes, p_e / c_e, color=blue_gradient[k], + label=f"C extension with output size: {size}", ) axs[0].grid() axs[1].grid() -axs[0].legend() +axs[1].legend(loc="upper left") axs[0].xaxis.set_major_formatter(ScalarFormatter(useMathText=False)) axs[1].xaxis.set_major_formatter(ScalarFormatter(useMathText=False)) @@ -84,8 +84,8 @@ axs[0].set_xlim(150, 4250) axs[1].set_xlim(150, 4250) -axs[0].set_ylim(0, 400) -axs[1].set_ylim(0, 400) +axs[0].set_ylim(0, 520) +axs[1].set_ylim(0, 520) axs[1].set_yticklabels([]) diff --git a/speed_test/sc_time_array.npy b/speed_test/sc_time_array.npy index 50f1d2d..70ccab6 100644 Binary files a/speed_test/sc_time_array.npy and b/speed_test/sc_time_array.npy differ diff --git a/speed_test/sc_time_array_with_err.npy b/speed_test/sc_time_array_with_err.npy index c8abfdc..d29de3c 100644 Binary files a/speed_test/sc_time_array_with_err.npy and b/speed_test/sc_time_array_with_err.npy differ diff --git a/speed_test/sn_time_array.npy b/speed_test/sn_time_array.npy index c1efe06..6db7389 100644 Binary files a/speed_test/sn_time_array.npy and b/speed_test/sn_time_array.npy differ diff --git a/speed_test/sn_time_array_with_err.npy b/speed_test/sn_time_array_with_err.npy index 0789c17..803d6a1 100644 Binary files a/speed_test/sn_time_array_with_err.npy and b/speed_test/sn_time_array_with_err.npy differ diff --git a/speed_test/sp_time_array.npy b/speed_test/sp_time_array.npy index 8182ab5..2dba487 100644 Binary files a/speed_test/sp_time_array.npy and b/speed_test/sp_time_array.npy differ diff --git a/speed_test/sp_time_array_with_err.npy b/speed_test/sp_time_array_with_err.npy index fb6728e..bf03da4 100644 Binary files a/speed_test/sp_time_array_with_err.npy and b/speed_test/sp_time_array_with_err.npy differ diff --git a/speed_test/speed_test.png b/speed_test/speed_test.png index b7c6d41..b3eb7ae 100644 Binary files a/speed_test/speed_test.png and b/speed_test/speed_test.png differ diff --git a/speed_test/speed_test.py b/speed_test/speed_test.py index 907fc83..e68d3ac 100644 --- a/speed_test/speed_test.py +++ b/speed_test/speed_test.py @@ -1,9 +1,8 @@ from timeit import Timer import numpy as np -from matplotlib import pyplot as plt -from spectres import spectres as sp -from spectres import spectres_numba as sn +from spectres.spectral_resampling import spectres as sp +from spectres.spectral_resampling_numba import spectres_numba as sn from spectresc import spectres as sc @@ -113,7 +112,7 @@ def call_spectresc_err(size_i, size_o): timer = Timer( stmt=f"call_spectresn({size_i}, {size_o})", setup=( - "from spectres import spectres_numba as sn; from __main__" + "from spectres.spectral_resampling_numba import spectres_numba as sn; from __main__" " import call_spectresn" ), ) @@ -193,7 +192,7 @@ def call_spectresc_err(size_i, size_o): timer = Timer( stmt=f"call_spectresn_err({size_i}, {size_o})", setup=( - "from spectres import spectres_numba as sn; from __main__" + "from spectres.spectral_resampling_numba import spectres_numba as sn; from __main__" " import call_spectresn_err" ), ) diff --git a/speed_test/speed_test_20230425.png b/speed_test/speed_test_20230425.png new file mode 100644 index 0000000..b7c6d41 Binary files /dev/null and b/speed_test/speed_test_20230425.png differ diff --git a/src/spectresc/__init__.py b/src/spectresc/__init__.py index 5e8aae6..b1358ed 100644 --- a/src/spectresc/__init__.py +++ b/src/spectresc/__init__.py @@ -1,3 +1,3 @@ -from .spectres import spectres +from .spectresc import spectres -__all__ = [spectres] +__all__ = ["spectres"] diff --git a/src/spectresc/spectres.c b/src/spectresc/spectres.c index 7ab6dd5..3825bab 100644 --- a/src/spectresc/spectres.c +++ b/src/spectresc/spectres.c @@ -6,10 +6,8 @@ #include #include -double *make_bins(double *wavs, int wavs_len) +void make_bins(double *wavs, int wavs_len, double *edges, double *widths) { - double *edges = (double *)malloc(sizeof(double) * (wavs_len + 1)); - edges[0] = wavs[0] - (wavs[1] - wavs[0]) / 2.0; edges[wavs_len] = wavs[wavs_len - 1] + (wavs[wavs_len - 1] - wavs[wavs_len - 2]) / 2.0; @@ -18,7 +16,11 @@ double *make_bins(double *wavs, int wavs_len) edges[i] = (wavs[i] + wavs[i - 1]) / 2.0; } - return edges; + for (int i = 0; i < wavs_len - 1; i++) + { + widths[i] = edges[i + 1] - edges[i]; + } + widths[wavs_len - 1] = edges[wavs_len] - edges[wavs_len - 1]; } // Define the spectres function @@ -49,196 +51,192 @@ static PyObject *spectres(PyObject *self, PyObject *args, PyObject *kwargs) return NULL; } - // Create new arrays if input arrays are not contiguous - if (!PyArray_ISCONTIGUOUS(new_wavs_array)) - { - new_wavs_array = (PyArrayObject *)PyArray_Cast(new_wavs_array, NPY_DOUBLE); - if (!new_wavs_array) - return NULL; - } - if (!PyArray_ISCONTIGUOUS(spec_wavs_array)) - { - spec_wavs_array = (PyArrayObject *)PyArray_Cast(spec_wavs_array, NPY_DOUBLE); - if (!spec_wavs_array) - return NULL; - } - if (!PyArray_ISCONTIGUOUS(spec_fluxes_array)) - { - spec_fluxes_array = (PyArrayObject *)PyArray_Cast(spec_fluxes_array, NPY_DOUBLE); - if (!spec_fluxes_array) - return NULL; - } - double *new_wavs = (double *)PyArray_DATA(new_wavs_array); double *spec_wavs = (double *)PyArray_DATA(spec_wavs_array); double *spec_fluxes = (double *)PyArray_DATA(spec_fluxes_array); - // Get the length of the input arrays int new_wavs_len = (int)PyArray_DIM(new_wavs_array, 0); int spec_wavs_len = (int)PyArray_DIM(spec_wavs_array, 0); - double *spec_edges = make_bins(spec_wavs, spec_wavs_len); - double *new_edges = make_bins(new_wavs, new_wavs_len); - double *spec_widths = (double *)malloc(sizeof(double) * spec_wavs_len); + /* --- Determine flux array dimensions --- */ + int ndim = PyArray_NDIM(spec_fluxes_array); + npy_intp *flux_shape = PyArray_DIMS(spec_fluxes_array); - for (int i = 0; i < spec_wavs_len; i++) + if (flux_shape[ndim - 1] != spec_wavs_len) { - spec_widths[i] = spec_edges[i + 1] - spec_edges[i]; + PyErr_SetString(PyExc_ValueError, "Last dimension of spec_fluxes must match length of spec_wavs."); + return NULL; } - // Create empty arrays for populating resampled flux and error - double *new_fluxes = (double *)malloc(sizeof(double) * new_wavs_len); + int num_spectra = 1; + for (int i = 0; i < ndim - 1; i++) + { + num_spectra *= flux_shape[i]; + } + /* --- Make bins --- */ + double *spec_edges = malloc((spec_wavs_len + 1) * sizeof(double)); + double *spec_widths = malloc(spec_wavs_len * sizeof(double)); + make_bins(spec_wavs, spec_wavs_len, spec_edges, spec_widths); + + double *new_edges = malloc((new_wavs_len + 1) * sizeof(double)); + double *new_widths = malloc(new_wavs_len * sizeof(double)); + make_bins(new_wavs, new_wavs_len, new_edges, new_widths); + + /* --- Handle optional errors --- */ PyArrayObject *spec_errs_array = NULL; double *spec_errs = NULL; - double *new_errs = NULL; - if (spec_errs_obj != NULL) + if (spec_errs_obj != NULL && spec_errs_obj != Py_None) { spec_errs_array = (PyArrayObject *)PyArray_FROM_OTF(spec_errs_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + if (!spec_errs_array) + { + Py_XDECREF(new_wavs_array); + Py_XDECREF(spec_wavs_array); + Py_XDECREF(spec_fluxes_array); + return NULL; + } - if (!spec_errs_array || !PyArray_ISCONTIGUOUS(spec_errs_array)) + npy_intp *err_shape = PyArray_DIMS(spec_errs_array); + for (int i = 0; i < ndim; i++) { - if (spec_errs_array) - spec_errs_array = (PyArrayObject *)PyArray_Cast(spec_errs_array, NPY_DOUBLE); - if (!spec_errs_array) + if (err_shape[i] != flux_shape[i]) { - Py_XDECREF(new_wavs_array); - Py_XDECREF(spec_wavs_array); - Py_XDECREF(spec_fluxes_array); + PyErr_SetString(PyExc_ValueError, "spec_errs must have the same shape as spec_fluxes."); return NULL; } } spec_errs = (double *)PyArray_DATA(spec_errs_array); - new_errs = (double *)malloc(sizeof(double) * new_wavs_len); } - int start = 0, stop = 0, warned = 0; - for (int i = 0; i < new_wavs_len; i++) + /* --- Create output arrays --- */ + npy_intp *new_flux_shape = malloc(ndim * sizeof(npy_intp)); + memcpy(new_flux_shape, flux_shape, ndim * sizeof(npy_intp)); + new_flux_shape[ndim - 1] = new_wavs_len; + + PyArrayObject *new_fluxes_array = (PyArrayObject *)PyArray_SimpleNew(ndim, new_flux_shape, NPY_DOUBLE); + double *new_fluxes_data = (double *)PyArray_DATA(new_fluxes_array); + + PyArrayObject *new_errs_array = NULL; + double *new_errs_data = NULL; + if (spec_errs != NULL) { - if (new_edges[i] < spec_edges[0] || new_edges[i + 1] > spec_edges[spec_wavs_len]) - { - new_fluxes[i] = fill; - if (spec_errs != NULL) - { - new_errs[i] = fill; - } - if ((i == 0 || i == new_wavs_len - 1) && verbose && !warned) - { - warned = 1; - printf("SpectResC: new_wavs contains values outside the range " - "in spec_wavs, new_fluxes and new_errs will be filled " - "with the value set in the 'fill' keyword argument.\n"); - } - } - else + new_errs_array = (PyArrayObject *)PyArray_SimpleNew(ndim, new_flux_shape, NPY_DOUBLE); + new_errs_data = (double *)PyArray_DATA(new_errs_array); + } + + free(new_flux_shape); + + /* --- Loop over each spectrum --- */ + for (int s = 0; s < num_spectra; s++) + { + double *this_flux = spec_fluxes + s * spec_wavs_len; + double *this_err = spec_errs != NULL ? spec_errs + s * spec_wavs_len : NULL; + double *out_flux = new_fluxes_data + s * new_wavs_len; + double *out_err = spec_errs != NULL ? new_errs_data + s * new_wavs_len : NULL; + + int start = 0, stop = 0, warned = 0; + + for (int i = 0; i < new_wavs_len; i++) { - while (spec_edges[start + 1] <= new_edges[i]) + if (new_edges[i] < spec_edges[0] || new_edges[i + 1] > spec_edges[spec_wavs_len]) { - start += 1; + out_flux[i] = fill; + if (out_err != NULL) + out_err[i] = fill; + + if ((i == 0 || i == new_wavs_len - 1) && verbose && !warned) + { + warned = 1; + PyErr_WarnEx(PyExc_RuntimeWarning, + "SpectResC: new_wavs contains values outside the range in spec_wavs; filled with 'fill'.", + 1); + } + continue; } + while (spec_edges[start + 1] <= new_edges[i]) + start++; while (spec_edges[stop + 1] < new_edges[i + 1]) - { - stop += 1; - } + stop++; if (stop == start) { - new_fluxes[i] = spec_fluxes[start]; - if (spec_errs != NULL) - { - new_errs[i] = spec_errs[start]; - } + out_flux[i] = this_flux[start]; + if (out_err != NULL) + out_err[i] = this_err[start]; } else { double start_factor = ((spec_edges[start + 1] - new_edges[i]) / (spec_edges[start + 1] - spec_edges[start])); double end_factor = ((new_edges[i + 1] - spec_edges[stop]) / (spec_edges[stop + 1] - spec_edges[stop])); - spec_widths[start] *= start_factor; - spec_widths[stop] *= end_factor; - - // Populate new_fluxes spectrum and uncertainty arrays double f_widths_sum = 0.0; double spec_widths_sum = 0.0; double e_wid_sum = 0.0; - for (int j = start; j <= stop; j++) - { - f_widths_sum += spec_widths[j] * spec_fluxes[j]; - spec_widths_sum += spec_widths[j]; - if (spec_errs != NULL) - { - e_wid_sum += pow(spec_widths[j] * spec_errs[j], 2); - } - } - new_fluxes[i] = f_widths_sum / spec_widths_sum; - - if (spec_errs != NULL) + for (int j = start; j <= stop; j++) { - new_errs[i] = sqrt(e_wid_sum) / spec_widths_sum; + double weight = spec_widths[j]; + if (j == start) + weight *= start_factor; + if (j == stop) + weight *= end_factor; + + f_widths_sum += weight * this_flux[j]; + spec_widths_sum += weight; + if (out_err != NULL) + e_wid_sum += pow(weight * this_err[j], 2); } - // Put back the old bin widths to their initial values - spec_widths[start] /= start_factor; - spec_widths[stop] /= end_factor; + out_flux[i] = f_widths_sum / spec_widths_sum; + if (out_err != NULL) + out_err[i] = sqrt(e_wid_sum) / spec_widths_sum; } } } - // Free the memory free(spec_edges); - free(new_edges); free(spec_widths); + free(new_edges); + free(new_widths); - // Create NumPy arrays to return the data to Python - npy_intp dims[1] = {new_wavs_len}; - PyObject *new_fluxes_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, new_fluxes); - PyArray_ENABLEFLAGS((PyArrayObject *)new_fluxes_array, NPY_ARRAY_OWNDATA); + Py_XDECREF(new_wavs_array); + Py_XDECREF(spec_wavs_array); + Py_XDECREF(spec_fluxes_array); + Py_XDECREF(spec_errs_array); if (spec_errs != NULL) { - PyObject *new_errs_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, new_errs); - PyArray_ENABLEFLAGS((PyArrayObject *)new_errs_array, NPY_ARRAY_OWNDATA); PyObject *result_list = PyList_New(0); - PyList_Append(result_list, new_fluxes_array); - PyList_Append(result_list, new_errs_array); - - Py_XDECREF(new_wavs_array); - Py_XDECREF(spec_wavs_array); - Py_XDECREF(spec_fluxes_array); - Py_XDECREF(spec_errs_array); - + PyList_Append(result_list, (PyObject *)new_fluxes_array); + PyList_Append(result_list, (PyObject *)new_errs_array); return result_list; } else { - Py_XDECREF(new_wavs_array); - Py_XDECREF(spec_wavs_array); - Py_XDECREF(spec_fluxes_array); - - return new_fluxes_array; + return (PyObject *)new_fluxes_array; } } // Define the module methods -static PyMethodDef SpectrescMethods[] = { +static PyMethodDef SpectresMethods[] = { {"spectres", (PyCFunction)spectres, METH_VARARGS | METH_KEYWORDS, "Resample a spectrum onto a new wavelength grid."}, {NULL, NULL, 0, NULL}}; // Define the module structure -static struct PyModuleDef spectrescmodule = { +static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, - "spectresc", // Module name + "spectresc", // Submodule name, it must match the Extension name in setup.py "Python extension module for the spectres function in C", // Module description -1, - SpectrescMethods}; + SpectresMethods}; // Define the module initialization function PyMODINIT_FUNC PyInit_spectresc(void) { import_array(); // Initialize NumPy - return PyModule_Create(&spectrescmodule); + return PyModule_Create(&moduledef); } diff --git a/test/test_with_numpy_arange.py b/test/test_with_numpy_arange.py index adc3be8..d248aad 100644 --- a/test/test_with_numpy_arange.py +++ b/test/test_with_numpy_arange.py @@ -3,6 +3,11 @@ from spectresc import spectres as sc +# =========================================================== +# One-dimensional flux arrays +# =========================================================== + + def test_with_random_arrays(): inwaves = np.linspace(4000.0, 10000.0, 1000) outwaves = np.linspace(4500.0, 8500.0, 731) @@ -80,3 +85,57 @@ def test_with_random_arrays_against_spectres_output_range_mismatched_right(): assert np.isclose(outfluxes_1_sp, outfluxes_1_sc, equal_nan=True).all() assert np.isclose(outfluxes_2_sp, outfluxes_2_sc, equal_nan=True).all() assert np.isclose(outfluxes_3_sp, outfluxes_3_sc, equal_nan=True).all() + + +# =========================================================== +# Multi-dimensional flux arrays +# =========================================================== + + +def test_multidim_fluxes_2d_against_spectres(): + """Test 2D flux array: shape (Nspec, Nwave)""" + inwaves = np.linspace(4000.0, 10000.0, 1000) + outwaves = np.linspace(4500.0, 8500.0, 731) + n_spec = 5 + influxes = np.random.random((n_spec, len(inwaves))) * 10000.0 + influxes_err = np.random.random((n_spec, len(inwaves))) * 100.0 + + outfluxes_sp, _ = sp(outwaves, inwaves, influxes, influxes_err) + outfluxes_sc, _ = sc(outwaves, inwaves, influxes, influxes_err) + + assert np.isclose(outfluxes_sp, outfluxes_sc).all() + + +def test_multidim_fluxes_3d_against_spectres(): + """Test 3D flux array: shape (Nbatch, Nspec, Nwave)""" + inwaves = np.linspace(4000.0, 10000.0, 1000) + outwaves = np.linspace(4500.0, 8500.0, 731) + n_batch, n_spec = 3, 4 + influxes = np.random.random((n_batch, n_spec, len(inwaves))) * 10000.0 + influxes_err = np.random.random((n_batch, n_spec, len(inwaves))) * 100.0 + + outfluxes_sp, _ = sp(outwaves, inwaves, influxes, influxes_err) + outfluxes_sc, _ = sc(outwaves, inwaves, influxes, influxes_err) + + assert np.isclose(outfluxes_sp, outfluxes_sc).all() + + +def test_multidim_fluxes_sum_conservation(): + """Check flux conservation across multi-spectra 2D array""" + inwaves = np.linspace(4000.0, 10000.0, 1000) + outwaves = np.linspace(4500.0, 8500.0, 731) + influxes = np.random.random((10, len(inwaves))) * 1000.0 + + outfluxes = sc(outwaves, inwaves, influxes) + waves_mask_for_inwaves = (inwaves >= outwaves[0]) & (inwaves <= outwaves[-1]) + + assert np.allclose( + np.sum(influxes[:, waves_mask_for_inwaves], axis=-1) * np.diff(inwaves).mean(), + np.sum(outfluxes, axis=-1) * np.diff(outwaves).mean(), + rtol=0.01, + ) + # the simple masking above will always give a larger sum on the outfluxes + assert ( + np.sum(influxes[:, waves_mask_for_inwaves], axis=-1) * np.diff(inwaves).mean() + < np.sum(outfluxes, axis=-1) * np.diff(outwaves).mean() + ).all()