diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 72f4042..c43f74a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -72,28 +72,20 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 + with: + fetch-depth: 0 # setuptools_scm needs history/tags for the version - name: Set up Python uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} - - name: Configure - run: cmake -S . -B build - - - name: Build - run: cmake --build build --parallel - - - name: Install Python package and dev dependencies - run: python -m pip install -e "python[dev]" + - name: Build and install the Python package (nanobind extension) + run: python -m pip install -e ".[dev]" - name: Run ruff if: ${{ matrix.python-version == '3.13' }} - working-directory: python run: python -m ruff check . - name: Run Python tests - working-directory: python - env: - LIBDEDX_SO: ${{ github.workspace }}/build/src/libdedx.so - run: python -m pytest tests -q + run: python -m pytest python/tests -q diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml new file mode 100644 index 0000000..7fa2756 --- /dev/null +++ b/.github/workflows/wheels.yml @@ -0,0 +1,101 @@ +name: Wheels + +on: + push: + branches: [main] + tags: ["v*"] + pull_request: + workflow_dispatch: + +concurrency: + group: wheels-${{ github.ref }} + cancel-in-progress: ${{ github.event_name == 'pull_request' }} + +# Least-privilege default for every job; publish jobs opt into id-token below. +permissions: + contents: read + +jobs: + build_wheels: + name: Wheels on ${{ matrix.os }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + # macOS arm64 + x86_64 are both built on the arm64 runner (see the + # [tool.cibuildwheel.macos] archs in pyproject.toml), so the scarce + # Intel macos-13 runner is not needed. + os: [ubuntu-latest, windows-latest, macos-14] + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 # needed for setuptools_scm to resolve the version + + - name: Build wheels + uses: pypa/cibuildwheel@v3.2.0 + + - uses: actions/upload-artifact@v4 + with: + name: cibw-wheels-${{ matrix.os }} + path: ./wheelhouse/*.whl + + build_sdist: + name: Build sdist + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 + + # `build` with no target builds the sdist and then a wheel *from that + # sdist*, which proves the sdist is self-contained. Only the sdist is + # published; the locally-built (unrepaired) wheel is discarded. + - name: Build sdist and verify it builds a wheel + run: pipx run build + + - uses: actions/upload-artifact@v4 + with: + name: cibw-sdist + path: dist/*.tar.gz + + # Dry-run publish to TestPyPI on every tag before the real PyPI release. + publish_testpypi: + name: Publish to TestPyPI + needs: [build_wheels, build_sdist] + runs-on: ubuntu-latest + if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v') + environment: + name: testpypi + url: https://test.pypi.org/p/libdedx + permissions: + id-token: write # OIDC trusted publishing + steps: + - uses: actions/download-artifact@v4 + with: + pattern: cibw-* + path: dist + merge-multiple: true + + - uses: pypa/gh-action-pypi-publish@release/v1 + with: + repository-url: https://test.pypi.org/legacy/ + skip-existing: true + + publish_pypi: + name: Publish to PyPI + needs: [publish_testpypi] + runs-on: ubuntu-latest + if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v') + environment: + name: pypi + url: https://pypi.org/p/libdedx + permissions: + id-token: write # OIDC trusted publishing + steps: + - uses: actions/download-artifact@v4 + with: + pattern: cibw-* + path: dist + merge-multiple: true + + - uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/.gitignore b/.gitignore index 46419a7..c1be312 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,12 @@ python/.pytest_cache/ python/*.egg-info/ *.pyc +# Python packaging (scikit-build-core / cibuildwheel) artifacts +dist/ +wheelhouse/ +*.egg-info/ +.pytest_cache/ + # Generated by CMake at configure time - do not commit libdedx/dedx_config.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 0063f62..20754e0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,6 +12,16 @@ include(CMakePackageConfigHelpers) option(DEDX_BUILD_EXAMPLES "Build libdedx example programs" ON) option(DEDX_BUILD_TESTS "Build libdedx test suite" ON) +option(DEDX_BUILD_PYTHON "Build the libdedx._core Python extension module" OFF) + +# When invoked by scikit-build-core (pip / cibuildwheel) build only the Python +# extension: skip the examples, tests, and the C-library install/packaging rules +# so the wheel stays a single self-contained module. +if(SKBUILD) + set(DEDX_BUILD_PYTHON ON) + set(DEDX_BUILD_EXAMPLES OFF) + set(DEDX_BUILD_TESTS OFF) +endif() # ---- Version from git tag ---- find_package(Git QUIET) @@ -24,6 +34,12 @@ if(GIT_FOUND) ERROR_QUIET ) endif() +# Under scikit-build-core prefer the version resolved by setuptools_scm so the C +# library and the Python package stay in lockstep (and so source builds without a +# .git directory still get a meaningful version). +if(SKBUILD AND DEFINED SKBUILD_PROJECT_VERSION_FULL AND NOT SKBUILD_PROJECT_VERSION_FULL STREQUAL "") + set(GIT_VERSION "${SKBUILD_PROJECT_VERSION_FULL}") +endif() if(NOT GIT_VERSION) set(GIT_VERSION "0.0.0-unknown") endif() @@ -50,7 +66,13 @@ endif() if(DEDX_BUILD_TESTS) add_subdirectory(tests) endif() +if(DEDX_BUILD_PYTHON) + add_subdirectory(python) +endif() +# The exported CMake package files and CPack packaging are only relevant for a +# regular C-library install, not for the Python wheel build. +if(NOT SKBUILD) configure_package_config_file( "${PROJECT_SOURCE_DIR}/cmake/dedxConfig.cmake.in" "${PROJECT_BINARY_DIR}/dedxConfig.cmake" @@ -126,3 +148,4 @@ configure_file( "${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake" IMMEDIATE @ONLY) add_custom_target(uninstall COMMAND ${CMAKE_COMMAND} -P ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake) +endif() # NOT SKBUILD diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 3d7d15d..24acd6d 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -67,3 +67,33 @@ locking. The intended fix is to audit the library for shared mutable state and either make workspaces single-thread-owned by design or add explicit synchronization where shared access is required. This is tracked as a known issue. + +The Python binding inherits this limitation: do not share a `Workspace` across +threads. (Calls do not release the GIL, so binding calls on a shared workspace +are serialized in practice, but the underlying library is still not thread-safe.) + +## Python binding + +The Python package lives under `python/` and is a +[nanobind](https://nanobind.readthedocs.io) extension (`libdedx._core`) built by +[scikit-build-core](https://scikit-build-core.readthedocs.io). It statically +links the `dedx` C target, so building it also compiles the C library. The whole +project is configured from the top-level `pyproject.toml`. + +```bash +pip install -e ".[dev]" # builds the extension in place +pytest python/tests +ruff check . +``` + +Notes for contributors: + +- The binding source (`python/src/dedx_core.cpp`) is C++17 and follows ordinary + C++ conventions; the C "declare variables at the top of the block" rule above + applies to the C library, not to this file. +- `dedx_config` owns the element arrays it is given and frees them in + `dedx_free_config()`, so any pointer handed to it from the binding must be + `malloc`'d. See the `Config` wrapper for how ownership and array lengths are + kept consistent. +- The package version comes from `setuptools_scm` (git tags) and is fed into the + C library so the two stay in lockstep. diff --git a/README.md b/README.md index 1047199..b1b80a3 100644 --- a/README.md +++ b/README.md @@ -64,6 +64,33 @@ dedx_free_workspace(ws, &err); See the [examples/](examples/) directory for more usage patterns. +## Python binding + +`libdedx` is also available as a Python package with the C library statically +linked in, so there is no separate install or build step and no runtime +dependency on a shared library: + +```bash +pip install libdedx +``` + +```python +import libdedx +from libdedx import _core as dedx + +stp = libdedx.simple_stp(dedx.PROTON, dedx.WATER, 100.0) # MeV cm² / g + +# Workspace/config object model for repeated evaluations +ws = dedx.Workspace() +cfg = dedx.Config() +cfg.program, cfg.ion, cfg.target = dedx.PSTAR, dedx.PROTON, dedx.WATER +ws.load(cfg) +stp = ws.stp(cfg, 100.0) +``` + +See [python/README.md](python/README.md) for the full Python API and the +development workflow. + ## Building Requires CMake 3.21+ and a C11 compiler. diff --git a/index.rst b/index.rst index c1936e7..6425760 100644 --- a/index.rst +++ b/index.rst @@ -189,6 +189,38 @@ Passing ``-1`` for the program, ion, or target slot lists the available values for that level. +*************** +Python binding +*************** + +libdedx ships a Python package built as a `nanobind +`_ extension with `scikit-build-core +`_. The C library is statically +linked into the extension and the stopping-power tables are embedded, so the +wheels are self-contained — ``pip install libdedx`` needs no local C build or +shared library. + +.. code-block:: python + + import libdedx + from libdedx import _core as dedx + + # one-shot lookup, mass stopping power in MeV cm^2 / g + stp = libdedx.simple_stp(dedx.PROTON, dedx.WATER, 100.0) + + # workspace/config object model for repeated evaluations + ws = dedx.Workspace() + cfg = dedx.Config() + cfg.program, cfg.ion, cfg.target = dedx.PSTAR, dedx.PROTON, dedx.WATER + ws.load(cfg) + stp = ws.stp(cfg, 100.0) + +The low-level ``libdedx._core`` module mirrors the C API: the workspace/config +object model, custom compounds, CSDA range, inverse stopping power / range, +unit conversion, composition and I-value accessors, and the program/ion/material +lists and names. The package version is kept in lockstep with the C library. + + ***** Notes ***** diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..4da5371 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,87 @@ +[build-system] +requires = [ + "scikit-build-core>0.10", + "nanobind>=2", + "setuptools-scm>=8", +] +build-backend = "scikit_build_core.build" + +[project] +name = "libdedx" +description = "Python binding for the libdedx charged-particle stopping-power library" +readme = "python/README.md" +requires-python = ">=3.9" +license = "GPL-3.0-or-later" +license-files = ["COPYING"] +authors = [ + { name = "Jakob Toftegaard" }, + { name = "Niels Bassler" }, + { name = "Leszek Grzanka" }, +] +dependencies = ["numpy"] +dynamic = ["version"] + +[project.urls] +Homepage = "https://github.com/APTG/libdedx" +Issues = "https://github.com/APTG/libdedx/issues" + +[project.optional-dependencies] +test = ["pytest>=8", "numpy"] +dev = ["ruff>=0.5.0", "pytest>=8", "numpy"] + +[tool.scikit-build] +minimum-version = "build-system.requires" +build-dir = "build/{wheel_tag}" +# The repo root is the CMake source tree; the extension is added when SKBUILD +# turns on DEDX_BUILD_PYTHON. Only the pure-Python package lives under python/. +wheel.packages = ["python/libdedx"] + +[tool.scikit-build.cmake] +version = ">=3.21" + +[tool.scikit-build.metadata.version] +provider = "scikit_build_core.metadata.setuptools_scm" + +[tool.scikit-build.sdist] +# Ship the C sources/headers needed to build the extension from an sdist. The +# stopping-power tables are already embedded in src/, so the raw data/ tree is +# not required to build a wheel. +include = [ + "CMakeLists.txt", + "include/**", + "src/**", + "python/**", +] + +[tool.setuptools_scm] +version_scheme = "post-release" +local_scheme = "node-and-date" + +[tool.pytest.ini_options] +testpaths = ["python/tests"] + +[tool.ruff] +line-length = 120 +target-version = "py39" +extend-exclude = ["build"] + +[tool.ruff.lint] +select = ["E", "F", "W"] + +[tool.cibuildwheel] +# CPython 3.9–3.14 on 64-bit platforms; PyPy and musllinux are skipped for now. +build = "cp39-* cp310-* cp311-* cp312-* cp313-* cp314-*" +skip = "pp* *-musllinux*" +build-frontend = "build" +# Install-and-import smoke test plus the full test suite for every wheel. +test-requires = "pytest numpy" +test-command = "pytest {project}/python/tests" + +[tool.cibuildwheel.macos] +# Build both architectures on a single (arm64) runner: arm64 natively and +# x86_64 by cross-compiling. This avoids GitHub's scarce Intel `macos-13` +# runners, which queue for a very long time. +archs = ["arm64", "x86_64"] + +[tool.cibuildwheel.macos.environment] +MACOSX_DEPLOYMENT_TARGET = "11.0" diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt new file mode 100644 index 0000000..138f6a5 --- /dev/null +++ b/python/CMakeLists.txt @@ -0,0 +1,43 @@ +# Builds the libdedx._core nanobind extension. This file is added by the +# top-level CMakeLists.txt when DEDX_BUILD_PYTHON is ON (scikit-build-core sets +# this automatically via SKBUILD). The extension statically links the `dedx` +# library target defined in src/CMakeLists.txt, producing a self-contained +# single-file module. + +enable_language(CXX) +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +# Locate the Python development headers and nanobind's CMake package. When built +# under scikit-build-core the interpreter and nanobind are provided by the build +# environment; the snippet below also works for a plain `cmake` invocation. +if(NOT Python_EXECUTABLE) + find_package(Python 3.9 REQUIRED COMPONENTS Interpreter Development.Module) +else() + find_package(Python 3.9 REQUIRED COMPONENTS Development.Module) +endif() + +execute_process( + COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir + OUTPUT_STRIP_TRAILING_WHITESPACE + OUTPUT_VARIABLE nanobind_ROOT +) +find_package(nanobind CONFIG REQUIRED) + +nanobind_add_module(_core + NB_STATIC + src/dedx_core.cpp +) +target_link_libraries(_core PRIVATE dedx) +# The binding calls a couple of internal helpers (e.g. dedx_internal_get_nucleon) +# whose headers live in src/, which the `dedx` target only exposes privately. +target_include_directories(_core PRIVATE "${PROJECT_SOURCE_DIR}/src") + +# Install the compiled module into the libdedx package directory inside the wheel. +# Python extension modules are MODULE libraries (LIBRARY artifact on every +# platform), but RUNTIME is listed too so the install is correct regardless of +# how a given generator/platform classifies the output. +install(TARGETS _core + LIBRARY DESTINATION libdedx + RUNTIME DESTINATION libdedx +) diff --git a/python/README.md b/python/README.md index 5bcb3fb..6bbb13f 100644 --- a/python/README.md +++ b/python/README.md @@ -1,7 +1,54 @@ # libdedx Python binding -This directory contains the `ctypes`-based Python binding for `libdedx`. +`libdedx` provides Python bindings for the +[libdedx](https://github.com/APTG/libdedx) charged-particle stopping-power +library. The bindings are a [nanobind](https://nanobind.readthedocs.io) +extension built with +[scikit-build-core](https://scikit-build-core.readthedocs.io). The libdedx C +library is **statically linked** into the extension and its data tables are +embedded, so wheels are self-contained — no separate shared library or system +install is required. -The Python module loads the native `libdedx` shared library at runtime. -Build and install the C library first, or set `LIBDEDX_SO` to the path of -the built shared library during development and testing. +## Installation + +```bash +pip install libdedx +``` + +## Quick start + +```python +import libdedx +from libdedx import _core as dedx + +# High-level convenience helpers +print(libdedx.get_version()) +stp = libdedx.get_stp(dedx.PSTAR, dedx.HYDROGEN, dedx.WATER_LIQUID, 100.0) # MeV cm2/g + +# Low-level workspace/config object model +ws = dedx.Workspace() +cfg = dedx.Config() +cfg.program = dedx.PSTAR +cfg.ion = dedx.HYDROGEN +cfg.target = dedx.WATER_LIQUID +ws.load(cfg) + +stp = ws.stp(cfg, 100.0) # mass stopping power, MeV cm2/g +rng = ws.csda(cfg, 100.0) # CSDA range, g/cm2 +energy = ws.inverse_csda(cfg, rng) # invert the range back to energy +``` + +Custom compounds, inverse stopping power, unit conversion, the +program/ion/material lists and names, composition and I-value accessors are all +exposed through `libdedx._core`. See `python/tests/` for further examples. + +## Development + +The whole project (C library + Python extension) is configured from the +top-level `pyproject.toml`. From the repository root: + +```bash +pip install -e ".[dev]" # builds the nanobind extension in place +pytest python/tests +ruff check . +``` diff --git a/python/libdedx/__init__.py b/python/libdedx/__init__.py index 2f7f17a..103116c 100644 --- a/python/libdedx/__init__.py +++ b/python/libdedx/__init__.py @@ -1,11 +1,106 @@ -"""Python package for the libdedx ctypes binding.""" +"""Python binding for the libdedx stopping-power library. -from ._api import get_csda_table, get_default_table, get_stp, get_stp_table, get_version +The compiled :mod:`libdedx._core` nanobind extension statically links the +libdedx C library, so importing this package requires no separately installed +shared library. It exposes both a faithful low-level API (the ``Workspace`` and +``Config`` object model plus the module-level functions) and a few high-level +convenience helpers that mirror the historical ctypes binding. +""" + +from __future__ import annotations + +import numpy as np + +from . import _core +from ._core import ( + Config, + Workspace, + composition, + convert_units, + csda_range_table, + default_energy_stp_table, + error_string, + i_value, + ion_list, + ion_name, + material_list, + material_name, + max_energy, + min_energy, + program_list, + program_name, + program_version, + simple_stp, + simple_stp_for_program, + stp_table, + stp_table_size, + version, + version_string, +) __all__ = [ + # low-level extension module + "_core", + # object model + "Config", + "Workspace", + # low-level functions + "composition", + "convert_units", + "csda_range_table", + "default_energy_stp_table", + "error_string", + "i_value", + "ion_list", + "ion_name", + "material_list", + "material_name", + "max_energy", + "min_energy", + "program_list", + "program_name", + "program_version", + "simple_stp", + "simple_stp_for_program", + "stp_table", + "stp_table_size", + "version", + "version_string", + # high-level convenience helpers "get_version", "get_stp", "get_stp_table", "get_default_table", "get_csda_table", ] + + +def get_version() -> str: + """Return the libdedx version as a ``major.minor.patch`` string.""" + major, minor, patch = version() + return f"{major}.{minor}.{patch}" + + +def get_stp(program: int, ion: int, target: int, energy: float) -> float: + """Return mass stopping power in MeV cm2/g for a single energy in MeV/nucl.""" + return simple_stp_for_program(program, ion, target, float(energy)) + + +def get_stp_table(program: int, ion: int, target: int, energies) -> np.ndarray: + """Return stopping powers (MeV cm2/g) for an array of energies (MeV/nucl).""" + energies = np.ascontiguousarray(energies, dtype=np.float64) + return stp_table(program, ion, target, energies) + + +def get_default_table(program: int, ion: int, target: int) -> tuple[np.ndarray, np.ndarray]: + """Return ``(energies, stps)`` for the built-in tabulated data points. + + ``energies`` are in MeV/nucl and ``stps`` in MeV cm2/g. + """ + return default_energy_stp_table(program, ion, target) + + +def get_csda_table(program: int, ion: int, target: int, energies) -> np.ndarray: + """Return CSDA ranges (g/cm2) for an array of energies (MeV/nucl).""" + energies = np.ascontiguousarray(energies, dtype=np.float64) + return csda_range_table(program, ion, target, energies) diff --git a/python/libdedx/__init__.pyi b/python/libdedx/__init__.pyi deleted file mode 100644 index 2635f7e..0000000 --- a/python/libdedx/__init__.pyi +++ /dev/null @@ -1,7 +0,0 @@ -def get_version() -> str: ... -def get_stp(program: int, ion: int, target: int, energy: float) -> float: ... -def get_stp_table(program: int, ion: int, target: int, energies: list[float]) -> list[float]: ... -def get_default_table(program: int, ion: int, target: int) -> tuple[list[float], list[float]]: ... -def get_csda_table(program: int, ion: int, target: int, energies: list[float]) -> list[float]: ... - -__all__: list[str] diff --git a/python/libdedx/_api.py b/python/libdedx/_api.py deleted file mode 100644 index 3dd9aa6..0000000 --- a/python/libdedx/_api.py +++ /dev/null @@ -1,120 +0,0 @@ -"""Public ctypes-backed API for libdedx.""" - -import ctypes -import ctypes.util -import os - - -def _load_library(): - path = os.environ.get("LIBDEDX_SO") - if path: - return ctypes.CDLL(path) - name = ctypes.util.find_library("dedx") - if name is None: - raise OSError( - "libdedx shared library not found. " - "Install libdedx or set LIBDEDX_SO to the path of the shared library." - ) - return ctypes.CDLL(name) - - -_lib = _load_library() - -_c_int_p = ctypes.POINTER(ctypes.c_int) -_c_float_p = ctypes.POINTER(ctypes.c_float) -_c_double_p = ctypes.POINTER(ctypes.c_double) - -_lib.dedx_get_simple_stp_for_program.argtypes = [ - ctypes.c_int, ctypes.c_int, ctypes.c_int, ctypes.c_float, _c_int_p, -] -_lib.dedx_get_simple_stp_for_program.restype = ctypes.c_float - -_lib.dedx_get_stp_table_size.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.c_int] -_lib.dedx_get_stp_table_size.restype = ctypes.c_int - -_lib.dedx_get_stp_table.argtypes = [ - ctypes.c_int, ctypes.c_int, ctypes.c_int, - ctypes.c_int, _c_float_p, _c_float_p, -] -_lib.dedx_get_stp_table.restype = ctypes.c_int - -_lib.dedx_fill_default_energy_stp_table.argtypes = [ - ctypes.c_int, ctypes.c_int, ctypes.c_int, _c_float_p, _c_float_p, -] -_lib.dedx_fill_default_energy_stp_table.restype = ctypes.c_int - -_lib.dedx_get_csda_range_table.argtypes = [ - ctypes.c_int, ctypes.c_int, ctypes.c_int, - ctypes.c_int, _c_float_p, _c_double_p, -] -_lib.dedx_get_csda_range_table.restype = ctypes.c_int - -_lib.dedx_get_version.argtypes = [_c_int_p, _c_int_p, _c_int_p] -_lib.dedx_get_version.restype = None - - -def _check(err): - if err.value != 0: - raise RuntimeError(f"libdedx error code {err.value}") - - -def get_version(): - """Return the libdedx version as a ``major.minor.patch`` string.""" - major = ctypes.c_int(0) - minor = ctypes.c_int(0) - patch = ctypes.c_int(0) - _lib.dedx_get_version(ctypes.byref(major), ctypes.byref(minor), ctypes.byref(patch)) - return f"{major.value}.{minor.value}.{patch.value}" - - -def get_stp(program, ion, target, energy): - """Return mass stopping power in MeV cm²/g for a single energy in MeV/nucl.""" - err = ctypes.c_int(0) - result = _lib.dedx_get_simple_stp_for_program( - program, ion, target, float(energy), ctypes.byref(err) - ) - _check(err) - return float(result) - - -def get_stp_table(program, ion, target, energies): - """Return stopping powers (MeV cm²/g) for a list of energies (MeV/nucl).""" - n = len(energies) - e_arr = (ctypes.c_float * n)(*energies) - s_arr = (ctypes.c_float * n)() - ret = _lib.dedx_get_stp_table(program, ion, target, n, e_arr, s_arr) - if ret != 0: - raise RuntimeError(f"libdedx error code {ret}") - return list(s_arr) - - -def get_default_table(program, ion, target): - """Return (energies, stps) for the built-in tabulated data points. - - energies -- list of energies in MeV/nucl - stps -- list of stopping powers in MeV cm^2/g - """ - n = _lib.dedx_get_stp_table_size(program, ion, target) - if n < 0: - raise RuntimeError(f"libdedx error code {n}") - if n == 0: - raise RuntimeError( - f"No tabulated data for program={program}, ion={ion}, target={target}" - ) - e_arr = (ctypes.c_float * n)() - s_arr = (ctypes.c_float * n)() - ret = _lib.dedx_fill_default_energy_stp_table(program, ion, target, e_arr, s_arr) - if ret < 0: - raise RuntimeError(f"libdedx error code {ret}") - return list(e_arr), list(s_arr) - - -def get_csda_table(program, ion, target, energies): - """Return CSDA ranges (g/cm^2) for a list of energies (MeV/nucl).""" - n = len(energies) - e_arr = (ctypes.c_float * n)(*energies) - r_arr = (ctypes.c_double * n)() - ret = _lib.dedx_get_csda_range_table(program, ion, target, n, e_arr, r_arr) - if ret != 0: - raise RuntimeError(f"libdedx error code {ret}") - return list(r_arr) diff --git a/python/libdedx/_core.pyi b/python/libdedx/_core.pyi new file mode 100644 index 0000000..20f3715 --- /dev/null +++ b/python/libdedx/_core.pyi @@ -0,0 +1,120 @@ +"""Type stubs for the libdedx._core nanobind extension.""" + +from typing import Optional + +import numpy as np +import numpy.typing as npt + +# --- program identifiers --- +ASTAR: int +PSTAR: int +ESTAR: int +MSTAR: int +ICRU73_OLD: int +ICRU73: int +ICRU49: int +ICRU: int +DEFAULT: int +BETHE_EXT00: int + +# --- aggregate states --- +DEFAULT_STATE: int +GAS: int +CONDENSED: int + +# --- MSTAR modes --- +MSTAR_MODE_A: int +MSTAR_MODE_B: int +MSTAR_MODE_G: int +MSTAR_MODE_H: int +MSTAR_MODE_C: int +MSTAR_MODE_D: int +MSTAR_MODE_DEFAULT: int + +# --- interpolation modes --- +INTERPOLATION_LOG_LOG: int +INTERPOLATION_LINEAR: int +INTERPOLATION_DEFAULT: int + +# --- stopping-power units --- +MEVCM2G: int +MEVCM: int +KEVUM: int + +# --- common ions / materials --- +HYDROGEN: int +PROTON: int +HELIUM: int +CARBON: int +WATER: int +WATER_LIQUID: int +WATER_VAPOR: int +AIR: int + +class Config: + def __init__(self) -> None: ... + program: int + target: int + ion: int + compound_state: int + interpolation_mode: int + mstar_mode: int + i_value: float + rho: float + elements_id: list[int] + elements_atoms: list[int] + elements_mass_fraction: list[float] + elements_i_value: list[float] + ion_a: int + @property + def cfg_id(self) -> int: ... + @property + def bragg_used(self) -> bool: ... + @property + def loaded(self) -> bool: ... + @property + def target_name(self) -> Optional[str]: ... + @property + def ion_name(self) -> Optional[str]: ... + @property + def program_name(self) -> Optional[str]: ... + +class Workspace: + def __init__(self, count: int = 1) -> None: ... + def load(self, config: Config) -> int: ... + def stp(self, config: Config, energy: float) -> float: ... + def csda(self, config: Config, energy: float) -> float: ... + def inverse_stp(self, config: Config, stp: float, side: int) -> float: ... + def inverse_csda(self, config: Config, range: float) -> float: ... + @property + def datasets(self) -> int: ... + @property + def active_datasets(self) -> int: ... + +def version() -> tuple[int, int, int]: ... +def version_string() -> str: ... +def error_string(err: int) -> str: ... +def program_name(program: int) -> str: ... +def program_version(program: int) -> str: ... +def material_name(material: int) -> str: ... +def ion_name(ion: int) -> str: ... +def program_list() -> list[int]: ... +def material_list(program: int) -> list[int]: ... +def ion_list(program: int) -> list[int]: ... +def min_energy(program: int, ion: int) -> float: ... +def max_energy(program: int, ion: int) -> float: ... +def i_value(target: int) -> float: ... +def composition(target: int) -> npt.NDArray[np.float64]: ... +def simple_stp(ion: int, target: int, energy: float) -> float: ... +def simple_stp_for_program(program: int, ion: int, target: int, energy: float) -> float: ... +def stp_table_size(program: int, ion: int, target: int) -> int: ... +def stp_table(program: int, ion: int, target: int, energies: npt.ArrayLike) -> npt.NDArray[np.float64]: ... +def csda_range_table( + program: int, ion: int, target: int, energies: npt.ArrayLike +) -> npt.NDArray[np.float64]: ... +def default_energy_stp_table( + program: int, ion: int, target: int +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: ... +def convert_units( + old_unit: int, new_unit: int, material: int, values: npt.ArrayLike +) -> npt.NDArray[np.float64]: ... diff --git a/python/pyproject.toml b/python/pyproject.toml deleted file mode 100644 index f8da0c0..0000000 --- a/python/pyproject.toml +++ /dev/null @@ -1,40 +0,0 @@ -[build-system] -requires = ["setuptools>=68", "wheel", "setuptools-scm>=8"] -build-backend = "setuptools.build_meta" - -[project] -name = "libdedx" -description = "Python ctypes binding for the libdedx stopping-power library" -readme = "README.md" -requires-python = ">=3.9" -license = {text = "LGPL-2.0-or-later"} -dynamic = ["version"] - -[project.optional-dependencies] -test = ["pytest>=8"] -dev = ["ruff>=0.5.0", "pytest>=8"] - -[tool.setuptools] -include-package-data = true - -[tool.setuptools.packages.find] -where = ["."] -include = ["libdedx"] - -[tool.setuptools.package-data] -libdedx = ["__init__.pyi", "py.typed"] - -[tool.pytest.ini_options] -testpaths = ["tests"] - -[tool.ruff] -line-length = 120 -target-version = "py39" - -[tool.ruff.lint] -select = ["E", "F", "W"] - -[tool.setuptools_scm] -root = ".." -version_scheme = "post-release" -local_scheme = "node-and-date" diff --git a/python/src/dedx_core.cpp b/python/src/dedx_core.cpp new file mode 100644 index 0000000..0187b64 --- /dev/null +++ b/python/src/dedx_core.cpp @@ -0,0 +1,532 @@ +// nanobind binding for libdedx — the low-level ``libdedx._core`` extension. +// +// This module statically links the libdedx C library and exposes the full +// public C API: the workspace/config object model, custom compounds, stopping +// power / CSDA range / inverse lookups, unit conversion, composition and +// I-value accessors, energy bounds, and the program/ion/material lists & names. +// +// numpy arrays are accepted (and returned) for the vectorised entry points so +// callers get array in / array out without manual buffer juggling. + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +extern "C" { +#include "dedx.h" +#include "dedx_periodic_table.h" +#include "dedx_tools.h" +#include "dedx_wrappers.h" +} + +namespace nb = nanobind; +using namespace nb::literals; + +namespace { + +// Raise a Python exception carrying the libdedx error string for a non-zero code. +void check(int err) { + if (err != 0) { + char buf[256]; + dedx_get_error_code(buf, err); + throw std::runtime_error("libdedx error " + std::to_string(err) + ": " + buf); + } +} + +// Read a -1 terminated static list returned by the C API into a vector. +std::vector read_terminated_list(const int *list) { + std::vector out; + if (list == nullptr) + return out; + for (int i = 0; list[i] != -1; ++i) + out.push_back(list[i]); + return out; +} + +// 1-D contiguous float64 input array (numpy's default dtype). +using FloatArrayIn = nb::ndarray, nb::c_contig, nb::device::cpu>; + +std::vector to_float_vector(const FloatArrayIn &arr) { + size_t n = arr.shape(0); + const double *src = arr.data(); + std::vector out(n); + for (size_t i = 0; i < n; ++i) + out[i] = static_cast(src[i]); + return out; +} + +// Build an owning 1-D numpy float64 array from a std::vector. +nb::ndarray make_numpy(std::vector &&values) { + auto *data = new std::vector(std::move(values)); + nb::capsule owner(data, [](void *p) noexcept { delete static_cast *>(p); }); + size_t n = data->size(); + return nb::ndarray(data->data(), {n}, owner); +} + +// --------------------------------------------------------------------------- +// Config: owns a heap dedx_config. The C library frees the element arrays and +// the struct itself via dedx_free_config(), so every pointer we install must be +// malloc()'d and ownership handed to the library. +// --------------------------------------------------------------------------- +class Config { + public: + Config() { + cfg_ = static_cast(std::calloc(1, sizeof(dedx_config))); + if (cfg_ == nullptr) + throw std::bad_alloc(); + cfg_->interpolation_mode = DEDX_INTERPOLATION_DEFAULT; + cfg_->mstar_mode = DEDX_MSTAR_MODE_DEFAULT; + } + + ~Config() { + int err = 0; + dedx_free_config(cfg_, &err); + } + + Config(const Config &) = delete; + Config &operator=(const Config &) = delete; + + dedx_config *raw() { return cfg_; } + + // --- scalar accessors ------------------------------------------------- + int program() const { return cfg_->program; } + void set_program(int v) { cfg_->program = v; } + int target() const { return cfg_->target; } + void set_target(int v) { cfg_->target = v; } + int ion() const { return cfg_->ion; } + void set_ion(int v) { cfg_->ion = v; } + int compound_state() const { return cfg_->compound_state; } + void set_compound_state(int v) { cfg_->compound_state = v; } + int interpolation_mode() const { return cfg_->interpolation_mode; } + void set_interpolation_mode(int v) { cfg_->interpolation_mode = v; } + int mstar_mode() const { return static_cast(cfg_->mstar_mode); } + void set_mstar_mode(int v) { cfg_->mstar_mode = static_cast(v); } + float i_value() const { return cfg_->i_value; } + void set_i_value(float v) { cfg_->i_value = v; } + float rho() const { return cfg_->rho; } + void set_rho(float v) { cfg_->rho = v; } + + // --- resolved fields ------------------------------------------------- + int cfg_id() const { return cfg_->cfg_id; } + int ion_a() const { return cfg_->ion_a; } + // Nucleon number; normally filled by load(), but writable to support + // isotopes or ions for which the default cannot be resolved. + void set_ion_a(int v) { cfg_->ion_a = v; } + bool bragg_used() const { return cfg_->bragg_used != 0; } + bool loaded() const { return cfg_->loaded != 0; } + + nb::object target_name() const { return name_or_none(cfg_->target_name); } + nb::object ion_name() const { return name_or_none(cfg_->ion_name); } + nb::object program_name() const { return name_or_none(cfg_->program_name); } + + // --- custom compound element arrays ---------------------------------- + std::vector elements_id() const { + return copy_ints(cfg_->elements_id, cfg_->elements_length); + } + void set_elements_id(const std::vector &v) { + // elements_id defines the compound length. If it changes, any previously + // set per-element arrays no longer match the new length and would cause + // out-of-bounds reads in dedx_load_config(), so drop them. + if (v.size() != cfg_->elements_length) + clear_dependent_arrays(); + replace_ints(&cfg_->elements_id, v); + cfg_->elements_length = static_cast(v.size()); + } + std::vector elements_atoms() const { + return copy_ints(cfg_->elements_atoms, cfg_->elements_length); + } + void set_elements_atoms(const std::vector &v) { + check_length(v.size()); + replace_ints(&cfg_->elements_atoms, v); + } + std::vector elements_mass_fraction() const { + return copy_floats(cfg_->elements_mass_fraction, cfg_->elements_length); + } + void set_elements_mass_fraction(const std::vector &v) { + check_length(v.size()); + replace_floats(&cfg_->elements_mass_fraction, v); + } + std::vector elements_i_value() const { + return copy_floats(cfg_->elements_i_value, cfg_->elements_length); + } + void set_elements_i_value(const std::vector &v) { + check_length(v.size()); + replace_floats(&cfg_->elements_i_value, v); + } + + private: + dedx_config *cfg_; + + static nb::object name_or_none(const char *s) { + if (s == nullptr) + return nb::none(); + return nb::cast(std::string(s)); + } + static std::vector copy_ints(const int *p, unsigned int n) { + std::vector out; + if (p != nullptr) + out.assign(p, p + n); + return out; + } + static std::vector copy_floats(const float *p, unsigned int n) { + std::vector out; + if (p != nullptr) + out.assign(p, p + n); + return out; + } + // Per-element arrays must be set *after* elements_id and match its length, + // otherwise the buffers handed to libdedx would be inconsistent with + // elements_length and read out of bounds. + void check_length(size_t n) const { + if (cfg_->elements_id == nullptr) + throw std::invalid_argument("set elements_id before the other element arrays"); + if (n != cfg_->elements_length) + throw std::invalid_argument("element array length must match elements_id length"); + } + void clear_dependent_arrays() { + std::free(cfg_->elements_atoms); + cfg_->elements_atoms = nullptr; + std::free(cfg_->elements_mass_fraction); + cfg_->elements_mass_fraction = nullptr; + std::free(cfg_->elements_i_value); + cfg_->elements_i_value = nullptr; + } + static void replace_ints(int **slot, const std::vector &v) { + std::free(*slot); + *slot = nullptr; + if (!v.empty()) { + *slot = static_cast(std::malloc(sizeof(int) * v.size())); + if (*slot == nullptr) + throw std::bad_alloc(); + std::memcpy(*slot, v.data(), sizeof(int) * v.size()); + } + } + static void replace_floats(float **slot, const std::vector &v) { + std::free(*slot); + *slot = nullptr; + if (!v.empty()) { + *slot = static_cast(std::malloc(sizeof(float) * v.size())); + if (*slot == nullptr) + throw std::bad_alloc(); + std::memcpy(*slot, v.data(), sizeof(float) * v.size()); + } + } +}; + +// --------------------------------------------------------------------------- +// Workspace: owns a dedx_workspace and drives evaluation against a Config. +// --------------------------------------------------------------------------- +class Workspace { + public: + explicit Workspace(unsigned int count) { + int err = 0; + ws_ = dedx_allocate_workspace(count, &err); + check(err); + if (ws_ == nullptr) + throw std::runtime_error("failed to allocate libdedx workspace"); + } + + ~Workspace() { + int err = 0; + dedx_free_workspace(ws_, &err); + } + + Workspace(const Workspace &) = delete; + Workspace &operator=(const Workspace &) = delete; + + int load(Config &config) { + int err = 0; + int id = dedx_load_config(ws_, config.raw(), &err); + check(err); + // load_config does not populate ion_a (only the C convenience wrappers + // do), yet csda()/inverse_*() require it. Fill it here so the object + // model "just works"; ignore failures for ions without nucleon data. + if (config.raw()->ion_a <= 0) { + int nerr = 0; + int a = dedx_internal_get_nucleon(config.raw()->ion, &nerr); + if (nerr == 0) + config.raw()->ion_a = a; + } + return id; + } + + float stp(Config &config, float energy) { + int err = 0; + float v = dedx_get_stp(ws_, config.raw(), energy, &err); + check(err); + return v; + } + + double csda(Config &config, float energy) { + int err = 0; + double v = dedx_get_csda(ws_, config.raw(), energy, &err); + check(err); + return v; + } + + double inverse_stp(Config &config, float stopping_power, int side) { + int err = 0; + double v = dedx_get_inverse_stp(ws_, config.raw(), stopping_power, side, &err); + check(err); + return v; + } + + double inverse_csda(Config &config, float range) { + int err = 0; + double v = dedx_get_inverse_csda(ws_, config.raw(), range, &err); + check(err); + return v; + } + + int datasets() const { return ws_->datasets; } + int active_datasets() const { return ws_->active_datasets; } + + private: + dedx_workspace *ws_; +}; + +// --------------------------------------------------------------------------- +// Free functions +// --------------------------------------------------------------------------- +std::tuple version() { + int major = 0, minor = 0, patch = 0; + dedx_get_version(&major, &minor, &patch); + return {major, minor, patch}; +} + +std::string version_string() { return dedx_get_version_string(); } + +std::string error_string(int err) { + char buf[256]; + dedx_get_error_code(buf, err); + return std::string(buf); +} + +std::string program_name(int program) { return dedx_get_program_name(program); } +std::string program_version(int program) { return dedx_get_program_version(program); } +std::string material_name(int material) { return dedx_get_material_name(material); } +std::string ion_name(int ion) { return dedx_get_ion_name(ion); } + +std::vector program_list() { return read_terminated_list(dedx_get_program_list()); } +std::vector material_list(int program) { return read_terminated_list(dedx_get_material_list(program)); } +std::vector ion_list(int program) { return read_terminated_list(dedx_get_ion_list(program)); } + +float min_energy(int program, int ion) { return dedx_get_min_energy(program, ion); } +float max_energy(int program, int ion) { return dedx_get_max_energy(program, ion); } + +float i_value(int target) { + int err = 0; + float v = dedx_get_i_value(target, &err); + check(err); + return v; +} + +// Returns an (N, 2) numpy array with [atomic_number, mass_fraction] per row. +nb::ndarray composition(int target) { + float buf[64][2]; + unsigned int len = 0; + int err = 0; + dedx_get_composition(target, buf, &len, &err); + check(err); + auto *data = new std::vector(static_cast(len) * 2); + for (unsigned int i = 0; i < len; ++i) { + (*data)[i * 2 + 0] = buf[i][0]; + (*data)[i * 2 + 1] = buf[i][1]; + } + nb::capsule owner(data, [](void *p) noexcept { delete static_cast *>(p); }); + return nb::ndarray(data->data(), {static_cast(len), 2}, owner); +} + +float simple_stp(int ion, int target, float energy) { + int err = 0; + float v = dedx_get_simple_stp(ion, target, energy, &err); + check(err); + return v; +} + +float simple_stp_for_program(int program, int ion, int target, float energy) { + int err = 0; + float v = dedx_get_simple_stp_for_program(program, ion, target, energy, &err); + check(err); + return v; +} + +int stp_table_size(int program, int ion, int target) { + return dedx_get_stp_table_size(program, ion, target); +} + +nb::ndarray stp_table(int program, int ion, int target, const FloatArrayIn &energies) { + std::vector e = to_float_vector(energies); + std::vector s(e.size()); + int ret = dedx_get_stp_table(program, ion, target, static_cast(e.size()), e.data(), s.data()); + check(ret); + std::vector out(s.begin(), s.end()); + return make_numpy(std::move(out)); +} + +nb::ndarray csda_range_table(int program, int ion, int target, const FloatArrayIn &energies) { + std::vector e = to_float_vector(energies); + std::vector r(e.size()); + int ret = dedx_get_csda_range_table(program, ion, target, static_cast(e.size()), e.data(), r.data()); + check(ret); + return make_numpy(std::move(r)); +} + +std::tuple, nb::ndarray> +default_energy_stp_table(int program, int ion, int target) { + int n = dedx_get_stp_table_size(program, ion, target); + if (n < 0) + check(n); + if (n == 0) + throw std::runtime_error("no tabulated data for the requested program/ion/target"); + std::vector e(n), s(n); + int ret = dedx_fill_default_energy_stp_table(program, ion, target, e.data(), s.data()); + if (ret < 0) + check(ret); + std::vector ed(e.begin(), e.end()); + std::vector sd(s.begin(), s.end()); + return {make_numpy(std::move(ed)), make_numpy(std::move(sd))}; +} + +nb::ndarray +convert_units_py(int old_unit, int new_unit, int material, const FloatArrayIn &values) { + std::vector in = to_float_vector(values); + // The C convert_units() short-circuits when the units match and leaves the + // output untouched, so handle that case here to return the values unchanged. + if (old_unit == new_unit) { + std::vector res(in.begin(), in.end()); + return make_numpy(std::move(res)); + } + std::vector out(in.size()); + int ret = convert_units(old_unit, new_unit, material, static_cast(in.size()), in.data(), out.data()); + check(ret); + std::vector res(out.begin(), out.end()); + return make_numpy(std::move(res)); +} + +} // namespace + +NB_MODULE(_core, m) { + m.doc() = "Low-level nanobind binding for the libdedx stopping-power C library."; + + // ---- enums / identifiers ------------------------------------------- + // Programs + m.attr("ASTAR") = (int)DEDX_ASTAR; + m.attr("PSTAR") = (int)DEDX_PSTAR; + m.attr("ESTAR") = (int)DEDX_ESTAR; + m.attr("MSTAR") = (int)DEDX_MSTAR; + m.attr("ICRU73_OLD") = (int)DEDX_ICRU73_OLD; + m.attr("ICRU73") = (int)DEDX_ICRU73; + m.attr("ICRU49") = (int)DEDX_ICRU49; + m.attr("ICRU") = (int)DEDX_ICRU; + m.attr("DEFAULT") = (int)DEDX_DEFAULT; + m.attr("BETHE_EXT00") = (int)DEDX_BETHE_EXT00; + + // Aggregate states + m.attr("DEFAULT_STATE") = (int)DEDX_DEFAULT_STATE; + m.attr("GAS") = (int)DEDX_GAS; + m.attr("CONDENSED") = (int)DEDX_CONDENSED; + + // MSTAR modes + m.attr("MSTAR_MODE_A") = (int)DEDX_MSTAR_MODE_A; + m.attr("MSTAR_MODE_B") = (int)DEDX_MSTAR_MODE_B; + m.attr("MSTAR_MODE_G") = (int)DEDX_MSTAR_MODE_G; + m.attr("MSTAR_MODE_H") = (int)DEDX_MSTAR_MODE_H; + m.attr("MSTAR_MODE_C") = (int)DEDX_MSTAR_MODE_C; + m.attr("MSTAR_MODE_D") = (int)DEDX_MSTAR_MODE_D; + m.attr("MSTAR_MODE_DEFAULT") = (int)DEDX_MSTAR_MODE_DEFAULT; + + // Interpolation modes + m.attr("INTERPOLATION_LOG_LOG") = (int)DEDX_INTERPOLATION_LOG_LOG; + m.attr("INTERPOLATION_LINEAR") = (int)DEDX_INTERPOLATION_LINEAR; + m.attr("INTERPOLATION_DEFAULT") = (int)DEDX_INTERPOLATION_DEFAULT; + + // Stopping-power units + m.attr("MEVCM2G") = (int)DEDX_MEVCM2G; + m.attr("MEVCM") = (int)DEDX_MEVCM; + m.attr("KEVUM") = (int)DEDX_KEVUM; + + // A few common ions / materials for convenience. + m.attr("HYDROGEN") = (int)DEDX_HYDROGEN; + m.attr("PROTON") = (int)DEDX_PROTON; + m.attr("HELIUM") = (int)DEDX_HELIUM; + m.attr("CARBON") = (int)DEDX_CARBON; + m.attr("WATER") = (int)DEDX_WATER; + m.attr("WATER_LIQUID") = (int)DEDX_WATER_LIQUID; + m.attr("WATER_VAPOR") = (int)DEDX_WATER_VAPOR; + m.attr("AIR") = (int)DEDX_AIR; + + // ---- Config -------------------------------------------------------- + nb::class_(m, "Config", "Stopping-power calculation configuration.") + .def(nb::init<>()) + .def_prop_rw("program", &Config::program, &Config::set_program) + .def_prop_rw("target", &Config::target, &Config::set_target) + .def_prop_rw("ion", &Config::ion, &Config::set_ion) + .def_prop_rw("compound_state", &Config::compound_state, &Config::set_compound_state) + .def_prop_rw("interpolation_mode", &Config::interpolation_mode, &Config::set_interpolation_mode) + .def_prop_rw("mstar_mode", &Config::mstar_mode, &Config::set_mstar_mode) + .def_prop_rw("i_value", &Config::i_value, &Config::set_i_value) + .def_prop_rw("rho", &Config::rho, &Config::set_rho) + .def_prop_rw("elements_id", &Config::elements_id, &Config::set_elements_id) + .def_prop_rw("elements_atoms", &Config::elements_atoms, &Config::set_elements_atoms) + .def_prop_rw("elements_mass_fraction", &Config::elements_mass_fraction, &Config::set_elements_mass_fraction) + .def_prop_rw("elements_i_value", &Config::elements_i_value, &Config::set_elements_i_value) + .def_prop_ro("cfg_id", &Config::cfg_id) + .def_prop_rw("ion_a", &Config::ion_a, &Config::set_ion_a) + .def_prop_ro("bragg_used", &Config::bragg_used) + .def_prop_ro("loaded", &Config::loaded) + .def_prop_ro("target_name", &Config::target_name) + .def_prop_ro("ion_name", &Config::ion_name) + .def_prop_ro("program_name", &Config::program_name); + + // ---- Workspace ----------------------------------------------------- + nb::class_(m, "Workspace", "Workspace holding preloaded stopping-power datasets.") + .def(nb::init(), "count"_a = 1) + .def("load", &Workspace::load, "config"_a, + "Load a configuration; returns the dataset id and populates resolved fields.") + .def("stp", &Workspace::stp, "config"_a, "energy"_a, + "Mass stopping power (MeV cm2/g) at the given energy (MeV/nucl).") + .def("csda", &Workspace::csda, "config"_a, "energy"_a, + "CSDA range (g/cm2) at the given energy (MeV/nucl). Requires a prior load().") + .def("inverse_stp", &Workspace::inverse_stp, "config"_a, "stp"_a, "side"_a, + "Energy (MeV/nucl) for a stopping power; side<0 low-energy branch, side>=0 high-energy branch.") + .def("inverse_csda", &Workspace::inverse_csda, "config"_a, "range"_a, + "Energy (MeV/nucl) for a CSDA range (g/cm2). Requires a prior load().") + .def_prop_ro("datasets", &Workspace::datasets) + .def_prop_ro("active_datasets", &Workspace::active_datasets); + + // ---- module-level functions ---------------------------------------- + m.def("version", &version, "Return the (major, minor, patch) library version."); + m.def("version_string", &version_string, "Return the full library version string."); + m.def("error_string", &error_string, "err"_a, "Human-readable description of an error code."); + m.def("program_name", &program_name, "program"_a); + m.def("program_version", &program_version, "program"_a); + m.def("material_name", &material_name, "material"_a); + m.def("ion_name", &ion_name, "ion"_a); + m.def("program_list", &program_list); + m.def("material_list", &material_list, "program"_a); + m.def("ion_list", &ion_list, "program"_a); + m.def("min_energy", &min_energy, "program"_a, "ion"_a); + m.def("max_energy", &max_energy, "program"_a, "ion"_a); + m.def("i_value", &i_value, "target"_a, "Mean excitation potential (eV) of a material."); + m.def("composition", &composition, "target"_a, + "Elemental composition as an (N, 2) array of [Z, mass_fraction]."); + m.def("simple_stp", &simple_stp, "ion"_a, "target"_a, "energy"_a); + m.def("simple_stp_for_program", &simple_stp_for_program, "program"_a, "ion"_a, "target"_a, "energy"_a); + m.def("stp_table_size", &stp_table_size, "program"_a, "ion"_a, "target"_a); + m.def("stp_table", &stp_table, "program"_a, "ion"_a, "target"_a, "energies"_a); + m.def("csda_range_table", &csda_range_table, "program"_a, "ion"_a, "target"_a, "energies"_a); + m.def("default_energy_stp_table", &default_energy_stp_table, "program"_a, "ion"_a, "target"_a, + "Return (energies, stps) for the built-in tabulated data points."); + m.def("convert_units", &convert_units_py, "old_unit"_a, "new_unit"_a, "material"_a, "values"_a, + "Convert an array of stopping-power values between unit systems."); +} diff --git a/python/tests/conftest.py b/python/tests/conftest.py deleted file mode 100644 index 3c91d9b..0000000 --- a/python/tests/conftest.py +++ /dev/null @@ -1,32 +0,0 @@ -import os -import sys -from pathlib import Path - - -def _candidate_library_paths(repo_root: Path) -> list[Path]: - if sys.platform.startswith("linux"): - lib_names = ["libdedx.so"] - elif sys.platform == "darwin": - lib_names = ["libdedx.dylib"] - elif os.name == "nt": - lib_names = ["dedx.dll", "libdedx.dll"] - else: - lib_names = ["libdedx.so"] - - candidates = [] - for build_dir in ("build", "build-debug", "build-release", "build-coverage"): - for lib_name in lib_names: - candidates.append(repo_root / build_dir / "src" / lib_name) - candidates.append(repo_root / build_dir / "libdedx" / lib_name) - return candidates - - -def pytest_configure() -> None: - if os.environ.get("LIBDEDX_SO"): - return - - repo_root = Path(__file__).resolve().parents[2] - for candidate in _candidate_library_paths(repo_root): - if candidate.exists(): - os.environ["LIBDEDX_SO"] = str(candidate) - return diff --git a/python/tests/test_core.py b/python/tests/test_core.py new file mode 100644 index 0000000..fa80d1b --- /dev/null +++ b/python/tests/test_core.py @@ -0,0 +1,163 @@ +"""Tests for the low-level libdedx._core API (workspace/config object model).""" + +import numpy as np +import pytest + +from libdedx import _core + +PSTAR = _core.PSTAR +HYDROGEN = _core.HYDROGEN +WATER = _core.WATER_LIQUID + + +def test_version_tuple(): + major, minor, patch = _core.version() + assert all(isinstance(v, int) for v in (major, minor, patch)) + + +def test_names_and_lists(): + assert isinstance(_core.program_name(PSTAR), str) + assert _core.program_name(PSTAR) + assert isinstance(_core.material_name(WATER), str) + assert isinstance(_core.ion_name(HYDROGEN), str) + + programs = _core.program_list() + assert PSTAR in programs + materials = _core.material_list(PSTAR) + assert len(materials) > 0 + ions = _core.ion_list(PSTAR) + assert HYDROGEN in ions + + +def test_energy_bounds(): + lo = _core.min_energy(PSTAR, HYDROGEN) + hi = _core.max_energy(PSTAR, HYDROGEN) + assert 0.0 < lo < hi + + +def test_accessors(): + iv = _core.i_value(WATER) + assert iv > 0.0 + comp = _core.composition(WATER) + assert comp.ndim == 2 and comp.shape[1] == 2 + assert comp.shape[0] >= 2 # water has at least H and O + # mass fractions should roughly sum to 1 + assert abs(comp[:, 1].sum() - 1.0) < 1e-2 + + +def test_workspace_config_roundtrip(): + ws = _core.Workspace(4) + cfg = _core.Config() + cfg.program = PSTAR + cfg.ion = HYDROGEN + cfg.target = WATER + ws.load(cfg) + assert cfg.loaded + assert cfg.cfg_id >= 0 + assert cfg.ion_a > 0 + assert cfg.program_name + assert cfg.target_name + + stp = ws.stp(cfg, 100.0) + assert stp > 0.0 + + +def test_csda_and_inverse(): + ws = _core.Workspace(8) + cfg = _core.Config() + cfg.program = PSTAR + cfg.ion = HYDROGEN + cfg.target = WATER + ws.load(cfg) + + energy = 100.0 + rng = ws.csda(cfg, energy) + assert rng > 0.0 + + recovered = ws.inverse_csda(cfg, rng) + assert recovered == pytest.approx(energy, rel=1e-2) + + stp = ws.stp(cfg, energy) + e_high = ws.inverse_stp(cfg, stp, side=1) + assert e_high == pytest.approx(energy, rel=5e-2) + + +def test_custom_compound_mass_fraction(): + # Water defined as a custom compound: H and O by mass fraction. + ws = _core.Workspace(2) + cfg = _core.Config() + cfg.program = PSTAR + cfg.ion = HYDROGEN + cfg.target = 0 # custom compound + cfg.elements_id = [1, 8] + cfg.elements_mass_fraction = [0.111894, 0.888106] + cfg.rho = 1.0 + ws.load(cfg) + assert cfg.loaded + assert cfg.bragg_used + + stp = ws.stp(cfg, 100.0) + assert stp > 0.0 + + +def test_custom_compound_by_atoms(): + # Water by atom counts (H2O). + ws = _core.Workspace(2) + cfg = _core.Config() + cfg.program = PSTAR + cfg.ion = HYDROGEN + cfg.target = 0 + cfg.elements_id = [1, 8] + cfg.elements_atoms = [2, 1] + cfg.rho = 1.0 + ws.load(cfg) + stp = ws.stp(cfg, 100.0) + assert stp > 0.0 + # mass fractions get derived from the atom counts during load + frac = cfg.elements_mass_fraction + assert len(frac) == 2 + assert frac[1] > frac[0] + + +def test_element_array_length_guards(): + cfg = _core.Config() + # Setting a per-element array before elements_id is rejected. + with pytest.raises(Exception): + cfg.elements_mass_fraction = [1.0] + + cfg.elements_id = [1, 8] + cfg.elements_mass_fraction = [0.11, 0.89] + # Mismatched length is rejected. + with pytest.raises(Exception): + cfg.elements_atoms = [2, 1, 3] + + # Shrinking elements_id drops the now-stale dependent arrays. + cfg.elements_id = [1] + assert cfg.elements_mass_fraction == [] + + +def test_convert_units(): + values = np.array([50.0, 25.0], dtype=np.float64) + # MeV cm2/g -> keV/um for liquid water (rho ~ 1 g/cm3): factor ~ 0.1 + converted = _core.convert_units(_core.MEVCM2G, _core.KEVUM, WATER, values) + assert converted.shape == values.shape + assert np.all(converted > 0.0) + # round trip + back = _core.convert_units(_core.KEVUM, _core.MEVCM2G, WATER, converted) + assert np.allclose(back, values, rtol=1e-4) + # identity conversion returns the values unchanged + same = _core.convert_units(_core.MEVCM2G, _core.MEVCM2G, WATER, values) + assert np.allclose(same, values) + + +def test_stp_table_matches_scalar(): + energies = np.array([1.0, 10.0, 100.0], dtype=np.float64) + table = _core.stp_table(PSTAR, HYDROGEN, WATER, energies) + assert table.shape == energies.shape + scalar = _core.simple_stp_for_program(PSTAR, HYDROGEN, WATER, 10.0) + assert table[1] == pytest.approx(scalar, rel=1e-5) + + +def test_error_string(): + assert isinstance(_core.error_string(0), str) + assert _core.error_string(201) # target not found diff --git a/python/tests/test_libdedx.py b/python/tests/test_libdedx.py index bb6399a..db5298d 100644 --- a/python/tests/test_libdedx.py +++ b/python/tests/test_libdedx.py @@ -1,23 +1,24 @@ -""" -Basic smoke tests for the libdedx Python binding. +"""High-level smoke tests for the libdedx Python binding. -Uses PSTAR (program=2), Hydrogen (ion=1), Liquid water (target=276) -as a reference combination with well-known stopping power values. +Uses PSTAR, hydrogen (Z=1) and liquid water as a reference combination with +well-known stopping power values. """ import libdedx -PROGRAM = 2 # PSTAR -ION = 1 # Hydrogen (Z=1) -TARGET = 276 # Liquid water +PROGRAM = libdedx._core.PSTAR +ION = libdedx._core.HYDROGEN +TARGET = libdedx._core.WATER_LIQUID -def test_get_version(capsys): +def test_get_version(): version = libdedx.get_version() - print(version) - captured = capsys.readouterr() assert version.count(".") == 2 - assert captured.out.strip() == version + + +def test_version_string_nonempty(): + assert isinstance(libdedx.version_string(), str) + assert libdedx.version_string() def test_get_stp_returns_positive(): @@ -45,8 +46,3 @@ def test_get_csda_table(): ranges = libdedx.get_csda_table(PROGRAM, ION, TARGET, energies) assert len(ranges) == len(energies) assert all(r > 0.0 for r in ranges) - - -# TODO: passing an invalid program number currently causes the library to -# segfault rather than returning an error code. A proper test can be added -# once input validation is implemented in the C library. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eb22244..665bcf8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,7 +42,9 @@ if(NOT WIN32) target_link_libraries(dedx PUBLIC m) endif() -# Shared library for Python ctypes and other dynamic consumers. +# Shared library for dynamic consumers. The Python wheel links the static `dedx` +# target instead, so the shared library is not built under SKBUILD. +if(NOT SKBUILD) add_library(dedx_shared SHARED $ ) @@ -68,7 +70,12 @@ endif() if(ANDROID) target_link_options(dedx_shared PRIVATE -Wl,-z,max-page-size=16384) endif() +endif() # NOT SKBUILD (dedx_shared) +# The C-library install rules are skipped for the Python wheel build (SKBUILD): +# the extension links `dedx` statically and the wheel should only contain the +# Python module, not the headers, archives, or exported CMake package. +if(NOT SKBUILD) install(FILES "${PROJECT_SOURCE_DIR}/include/dedx.h" "${PROJECT_SOURCE_DIR}/include/dedx_elements.h" @@ -90,3 +97,4 @@ install(EXPORT dedxTargets NAMESPACE dedx:: DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/dedx" ) +endif() # NOT SKBUILD