diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index a90ac36d..99c9ea1d 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -92,6 +92,9 @@ jobs: # dependency here. CIBW_ARCHS_LINUX: x86_64 CIBW_ARCHS_WINDOWS: AMD64 + # libgomp.so.1 is the one torch carries in torch/lib, so the wheel + # vendors none and the process has one OpenMP runtime. + CIBW_REPAIR_WHEEL_COMMAND_LINUX: auditwheel repair --exclude libgomp.so.1 -w {dest_dir} {wheel} # The library imports the system's DLLs and the OpenMP runtime torch # carries, so there is nothing to vendor and no repair step. CIBW_ENVIRONMENT_WINDOWS: CC=clang CXX=clang++ CMAKE_GENERATOR=Ninja @@ -125,7 +128,7 @@ jobs: grep -q 'licenses/external/finufft/LICENSE' contents.txt grep -q 'licenses/finufft-dependencies/DUCC0/LICENSE' contents.txt grep -q 'licenses/finufft-dependencies/xsimd/LICENSE' contents.txt - ! grep -Ei 'finufft/.*\.(so|dylib|dll)|libomp|libiomp' contents.txt + ! grep -Ei 'finufft/.*\.(so|dylib|dll)|libomp|libiomp|libgomp' contents.txt mkdir unpacked && python -m zipfile -e wheelhouse/*.whl unpacked case "${{ runner.os }}" in Linux) @@ -230,8 +233,9 @@ jobs: python -m build --wheel --outdir dist/ \ -C cmake.define.BARTORCH_CUDA=ON \ -C "cmake.define.BARTORCH_MATHDX_DIR=$MATHDX_DIR" - # The CUDA runtimes are excluded rather than vendored: torch already - # brings them, and bundling them would add hundreds of megabytes. + # The CUDA runtimes and libgomp are excluded rather than vendored: torch + # already brings them, bundling the first would add hundreds of + # megabytes, and a second libgomp would be a second OpenMP runtime. - name: Repair to a manylinux tag run: | /opt/python/cp312-cp312/bin/python -m pip install auditwheel @@ -239,14 +243,15 @@ jobs: /opt/python/cp312-cp312/bin/python -m auditwheel repair dist/*.whl \ -w wheelhouse/ --plat manylinux_2_28_x86_64 \ --exclude libcudart.so.12 --exclude libcufft.so.11 \ - --exclude libcublas.so.12 --exclude libcublasLt.so.12 + --exclude libcublas.so.12 --exclude libcublasLt.so.12 \ + --exclude libgomp.so.1 # cuFINUFFT is in the library, bound to the same cudart and cuFFT as - # BART's kernels, and nothing is vendored beside it but libgomp. + # BART's kernels, and nothing is vendored beside it. - name: Inspect the CUDA wheel run: | /opt/python/cp312-cp312/bin/python -m auditwheel show wheelhouse/*.whl /opt/python/cp312-cp312/bin/python -m zipfile -l wheelhouse/*.whl | tee contents.txt - ! grep -Ei 'libcudart|libcufft|libcublas|cufinufft.*\.so' contents.txt + ! grep -Ei 'libcudart|libcufft|libcublas|libgomp|cufinufft.*\.so' contents.txt grep -q 'licenses/finufft-dependencies/DUCC0/LICENSE' contents.txt # The wheel carries no interpreter, so it is loaded from the oldest and # the newest Python the package supports. diff --git a/AGENTS.md b/AGENTS.md index e7860e6c..c7b52431 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -222,18 +222,23 @@ as for BART's kernels. cuFINUFFT at this pin does not compile for sm_90 with CU takes for that architecture arrives later), so a CUDA build that includes 90 needs 12.1 or newer; the wheel is built with 12.8. -**One OpenMP runtime, and it is torch's where that matters.** +**One OpenMP runtime, and it is torch's.** `cmake/openmp.cmake` defines `OpenMP::OpenMP_C` and `OpenMP::OpenMP_CXX` for the whole build, BART linking the first and FINUFFT the second, and answers FINUFFT's own `find_package(OpenMP)` through `CMAKE_FIND_PACKAGE_REDIRECTS_DIR` -so it cannot pick a runtime of its own. On Linux it is the toolchain's -runtime. On Windows and macOS LLVM's runtime ends the process when a second +so it cannot pick a runtime of its own. On Linux the build links the +toolchain's libgomp, whose `libgomp.so.1` is also the name of the copy torch +carries in `torch/lib` from 2.7.1 on -- before that it had a hashed name, +hence the Linux floor in `pyproject.toml` -- so the copy torch has loaded +satisfies the library's `NEEDED` entry; the wheel is repaired with +`--exclude libgomp.so.1` and carries none, and `$ORIGIN/../torch/lib` is on +its runpath. On Windows and macOS LLVM's runtime ends the process when a second copy initialises (`OMP: Error #15`), and torch has already loaded one, so the library is linked against a stub that names torch's and lists the entry points clang calls: an import library for `libiomp5md.dll` made from `src/csrc/compat/libiomp5md.def`, and a text stub for `@rpath/libomp.dylib` written from the same list, with `@loader_path/../torch/lib` on the rpath. -`_lib._load` imports torch before it loads the library on both. clang 19 +`_lib._load` imports torch before it loads the library on every platform. clang 19 and later end a dynamically scheduled loop with `__kmpc_dispatch_deinit`, which the libomp in torch 2.3 to 2.5 does not export, so `src/csrc/substitute/openmp.c` answers it inside the library and forwards to @@ -242,7 +247,7 @@ macOS, hence the platform floor in `pyproject.toml`. macOS needs `omp.h` to compile, which Homebrew's `libomp` provides; nothing links against that copy. `tests/test_openmp.py` runs threaded torch and threaded FINUFFT in one process in either import order and reads the loaded images: one runtime -on macOS and Windows, and it is torch's, with no installed file changed. +on every platform, and it is torch's, with no installed file changed. `KMP_DUPLICATE_LIB_OK` is not used. **Underneath BART's own tools the seam is `nufft_create`, not the gridder.** @@ -314,7 +319,7 @@ BART's gridder along with everything else in that file. So `compute_psf`, `compute_psf2` and `compute_psf2_decomposed` are renamed too and `src/csrc/substitute/psf.c` answers to them: the same squared weights and basis, the same doubled grid, shifts and decomposition, with the transform in the middle being whichever -`nufft_create2` answers. `nlinv`, `moba`, `rtnlinv`, `noir/model2` and the +`nufft_create2` answers. `nlinv`, `moba`, `noir/model2` and the `psf` tool call these directly. What that is worth is in the numbers. On the un-doubled grid of a 16x16 @@ -612,8 +617,8 @@ points that read an array element by element rather than through `md_` -- taking a median of k-space -- are answered in `src/csrc/abi/host_reads.c` over a host copy of that one array, which is what BART already does for its own virtual pointers. What cannot be reached that way is a tool that allocates a temporary -of its own on the host and mixes it with its input: `pocsense` takes its -pattern from `md_alloc`, `nlinv` from `anon_cfl`, `ecalib` sets `bart_use_gpu` +of its own on the host and mixes it with its input: `nlinv` takes its +pattern from `anon_cfl`, `ecalib` sets `bart_use_gpu` from its own flag. `md_` operations take the host path unless every argument is on a device, and take it silently, so that is a segmentation fault rather than a slower answer. @@ -893,8 +898,8 @@ declared it. An input file that is not there is the one thing the command line names itself. A BART command that fails while loading its arguments leaves the library unable -to serve the next call in the same process -- `ecalib`, `nufft` and `pocsense` -handed a name with no file behind it all spin the call after them, while `fft` +to serve the next call in the same process -- `ecalib` and `nufft` +handed a name with no file behind it both spin the call after them, while `fft` does not -- so a caller who runs `main` twice would hang rather than see the second answer. `cli._missing` checks the names against the filesystem before BART is asked. @@ -926,10 +931,10 @@ reads no residual, and `iter2_pocs` is handed an `xupdate_op` that `pocs` never calls. What the sets are belongs to the projections, so `apps.pocsense` builds the application's three out of `linop.Sampling`, the range of a `linop.CartesianSense`, and a `priors` term conjugated by the transform -between the samples and the coil images. The three are bit-identical to -`tools.pocsense` on a grid, in two dimensions and in three, on an even grid -and on an odd one. The odd grid is what decides how that projection is -written: there the modulation is a phase rather than a sign, so the scaling +between the samples and the coil images, and `tests/test_apps.py` holds a +fully sampled phantom to be a fixed point of them and an answer to lie in the +range of the coils. The odd grid is what decides how the sparsity projection +is written: there the modulation is a phase rather than a sign, so the scaling has to ride in the same array BART puts it in rather than in a second multiply, and the multiply itself has to be BART's `md_zmul2` and `md_zmulc2` -- a `linop.Diagonal` and its adjoint -- rather than torch's. A complex @@ -1101,13 +1106,13 @@ reconstruction is held to its data rather than to its transform. Two exceptions: a caller who names an upsampling gets it -- `nufft_conf_s` carries two as BART's own default, so zero is what says nobody asked, and BART gets its two back before it sees the conf -- and the tools that calibrate before a -reconstruction is attempted (`nlinv`, `rtnlinv`, `ncalib`, in `_CALIBRATES`) +reconstruction is attempted (`nlinv` and `ncalib`, in `_CALIBRATES`) run at FINUFFT's own tolerance on the textbook grid, which costs a few megabytes at the resolution they fit sensitivities at. Every BART entry point that builds a NUFFT is served: `nufft` forward, -adjoint, inverse and Toeplitz, `pics` with and without a pattern, `sqpics`, -`nlinv`, `rtnlinv`, `moba`, `ncalib` and `linop.NUFFT`, on the host +adjoint, inverse and Toeplitz, `pics` with and without a pattern, +`nlinv`, `moba`, `ncalib` and `linop.NUFFT`, on the host and on the card, with BART's own gridder built zero times. `nlinv` and the network models build theirs against dimensions alone and hand the trajectory over afterwards, which is why `nufft_update_traj` installs one rather than @@ -1121,7 +1126,7 @@ index sets that are not axes (channels, parameter maps) are tuples too; no public argument takes a bitmask or a `-R` string. Hand-written wrappers convert their own; derived wrappers and what a curated one passes through by name follow `_call.TRANSLATED`, and `test_tools` fails on a dimension-like -argument that is in neither. `pics`, `sqpics`, `wshfl` and `moba` take +argument that is in neither. `pics`, `wshfl` and `moba` take `bartorch.priors` terms, serialized by `Regularizer._argument`. A flag's value can be an array rather than a number -- `pics(kspace, maps, t=traj)`, `-p` for a sampling pattern, `-B` for a basis -- and is registered and copied like any diff --git a/CMakeLists.txt b/CMakeLists.txt index 7533f908..741787ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,13 +26,14 @@ option(BARTORCH_CUDA "Build BART's CUDA kernels" OFF) set(BARTORCH_CUDA_ARCHITECTURES "75;80;86;89;90" CACHE STRING "Architectures to emit SASS for; PTX for the last one is added for newer cards") -set(CMAKE_C_STANDARD 11) +set(CMAKE_C_STANDARD 23) set(CMAKE_C_EXTENSIONS ON) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_EXTENSIONS OFF) set(CMAKE_POSITION_INDEPENDENT_CODE ON) set(CMAKE_C_VISIBILITY_PRESET hidden) set(CMAKE_CXX_VISIBILITY_PRESET hidden) +set(CMAKE_CUDA_VISIBILITY_PRESET hidden) set(CMAKE_VISIBILITY_INLINES_HIDDEN ON) if(NOT CMAKE_BUILD_TYPE) @@ -171,7 +172,7 @@ if(BARTORCH_CUDA) enable_language(CUDA) find_package(CUDAToolkit REQUIRED) - list(APPEND BART_C_DEFS USE_CUDA CUDA_GET_CUDA_DEVICE_NUM) + list(APPEND BART_C_DEFS USE_CUDA USE_GPU CUDA_GET_CUDA_DEVICE_NUM) set(CMAKE_CUDA_STANDARD 17) set(CMAKE_CUDA_STANDARD_REQUIRED ON) @@ -269,7 +270,7 @@ file(GLOB _tool_srcs CONFIGURE_DEPENDS "${BART_SRC}/*.c") list(APPEND BART_SOURCES ${_tool_srcs}) if(BARTORCH_CUDA) - foreach(_dir calib motion noncart num wavelet) + foreach(_dir calib motion noncart num simu wavelet) file(GLOB _cu CONFIGURE_DEPENDS "${BART_SRC}/${_dir}/*.cu") list(APPEND BART_SOURCES ${_cu}) endforeach() @@ -467,6 +468,10 @@ target_compile_definitions(bartorch PRIVATE # functions linked into this library rather than imported from its DLL # (misc/dllspec.h); and no named pipes, which are a POSIX file type. $<$:_FILE_OFFSET_BITS=64;BARTLIB_STATIC;NO_FIFO> + # Elsewhere misc/dllspec.h gives BARTLIB_API default visibility, which + # would put BART's seq_*, num_rand_init and debug_printf in the dynamic + # table beside the bartorch_* ABI; defining its guard leaves them hidden. + $<$>:BART_MISC_DLLSPEC;BARTLIB_API=;BARTLIB_CALL=> BARTORCH_BUILD_INFO="bart=${BART_VERSION},compiler=${CMAKE_C_COMPILER_ID}-${CMAKE_C_COMPILER_VERSION},nested=${BARTORCH_NESTED},openmp=${OpenMP_C_FOUND},cuda=${BARTORCH_CUDA},finufft=${BARTORCH_FINUFFT_VERSION},finufft_fft=${BARTORCH_FINUFFT_FFT_NAME},finufft_simd=${_finufft_simd_info},finufft_simd_mkl=${_finufft_mkl_info},long=${BARTORCH_SIZEOF_LONG}") target_compile_options(bartorch PRIVATE @@ -506,6 +511,13 @@ if(TARGET bartorch_mkl_fftw) set_property(TARGET bartorch APPEND PROPERTY INSTALL_RPATH "${BARTORCH_MKL_LIBRARY_DIR}") endif() +# Linux: libgomp.so.1 is torch's in a wheel, which carries none of its own; +# torch is its sibling package. The copy torch has loaded satisfies the name +# first, so this is where it is found when torch has not been imported. +if(NOT APPLE AND NOT WIN32 AND OpenMP_C_FOUND) + set_property(TARGET bartorch APPEND PROPERTY INSTALL_RPATH "$ORIGIN/../torch/lib") +endif() + # macOS: the OpenMP stub names @rpath/libomp.dylib, which is torch's. In a # wheel torch is the library's sibling package; an in-place build finds it # where the interpreter that configured the build has it. diff --git a/cmake/openmp.cmake b/cmake/openmp.cmake index c8eba376..1149e7d0 100644 --- a/cmake/openmp.cmake +++ b/cmake/openmp.cmake @@ -6,7 +6,9 @@ # of the library are bound to the same runtime. # # Linux: the toolchain's own runtime, found by CMake's FindOpenMP -- libgomp -# under GCC, libomp under clang. +# under GCC, libomp under clang. libgomp.so.1 is also the name of the copy +# torch carries, so the copy torch has loaded satisfies the library's NEEDED +# entry; the wheel is repaired with --exclude libgomp.so.1 and carries none. # # Windows and macOS: the runtime torch carries and has already loaded. LLVM's # runtime ends the process when a second copy of itself initialises (OMP: Error diff --git a/docs/api/apps.md b/docs/api/apps.md index 6afaf6cb..36261f14 100644 --- a/docs/api/apps.md +++ b/docs/api/apps.md @@ -17,5 +17,5 @@ runs an app in place of the BART command where one exists. | Object | BART command | Relation to the command | | --- | --- | --- | | {obj}`~bartorch.apps.pics` | `pics` | Identical output on a Cartesian grid; agrees to floating-point round-off along a trajectory | -| {obj}`~bartorch.apps.pocsense` | `pocsense` | Identical output | +| {obj}`~bartorch.apps.pocsense` | none | The POCSENSE projections swept by BART's `pocs` iteration | | {obj}`~bartorch.apps.mobafit` | `mobafit` | Same Gauss-Newton method over a TorchSim model; returns named parameter maps in physical units | diff --git a/docs/api/tools.md b/docs/api/tools.md index fbba7f60..14a27879 100644 --- a/docs/api/tools.md +++ b/docs/api/tools.md @@ -67,13 +67,10 @@ interface. | --- | --- | | {obj}`~bartorch.tools.pics` | Parallel-imaging compressed-sensing reconstruction (wrapped) | | {obj}`~bartorch.tools.nlinv` | Nonlinear inversion: image and sensitivities jointly (wrapped) | -| {obj}`~bartorch.tools.rtnlinv` | Real-time nonlinear inversion of a time series | | {obj}`~bartorch.tools.moba` | Model-based nonlinear reconstruction of parameter maps | | {obj}`~bartorch.tools.mobafit` | Voxel-wise fit of a signal model to contrast images | | {obj}`~bartorch.tools.looklocker` | $T_1$ from Look-Locker parameters $M_0$, $M_{ss}$ and $R_1^*$ | | {obj}`~bartorch.tools.itsense` | Iterative SENSE with $\ell_2$ regularization | -| {obj}`~bartorch.tools.sqpics` | Parallel-imaging compressed-sensing reconstruction, BART's `sqpics` variant of `pics` | -| {obj}`~bartorch.tools.pocsense` | POCSENSE reconstruction | | {obj}`~bartorch.tools.sake` | SAKE: low-rank matrix completion of k-space | | {obj}`~bartorch.tools.lrmatrix` | Multi-scale low-rank matrix completion | | {obj}`~bartorch.tools.homodyne` | Homodyne partial-Fourier reconstruction | diff --git a/docs/design/finufft-embedding.md b/docs/design/finufft-embedding.md index 983cb29c..04f80698 100644 --- a/docs/design/finufft-embedding.md +++ b/docs/design/finufft-embedding.md @@ -21,7 +21,7 @@ libbartorch (shared; exports the bartorch_* C ABI only) │ └── CUDA::cudart, CUDA::cufft ── dynamic, from the nvidia wheels torch brings └── OpenMP::OpenMP_C ────────────┤ └── cmake/openmp.cmake - Linux the toolchain's runtime (libgomp with GCC) + Linux libgomp.so.1 = torch/lib/libgomp.so.1 macOS @rpath/libomp.dylib = torch/lib/libomp.dylib Windows libiomp5md.dll, as torch loads it ``` @@ -110,9 +110,11 @@ runtime's own when the process has one. torch 2.2 carries no libomp on macOS, so the macOS floor is 2.3. `tests/test_openmp.py` is the check: threaded torch and threaded FINUFFT in -one process, both import orders, the loaded OpenMP images listed. On Linux -the one image found is torch's `libgomp.so.1`, which satisfies the library's -`NEEDED` entry as well. +one process, both import orders, the loaded OpenMP images listed, and one +image found on every platform. On Linux it is torch's `libgomp.so.1`, which +satisfies the library's `NEEDED` entry; the wheel is repaired with +`--exclude libgomp.so.1`, and torch's copy has carried that name from 2.7.1 +on, which is the Linux floor. ## FFT inside FINUFFT: DUCC0 and oneMKL diff --git a/docs/explanation/non-cartesian.md b/docs/explanation/non-cartesian.md index b611ab57..16997988 100644 --- a/docs/explanation/non-cartesian.md +++ b/docs/explanation/non-cartesian.md @@ -52,8 +52,8 @@ chooses its kernel width. | Setting | Default | Set by | | --- | --- | --- | -| Tolerance $\varepsilon$ | $10^{-3}$; $10^{-6}$ for `ncalib`, `nlinv` and `rtnlinv` | The kernel width: `width` of {class}`~bartorch.linop.NUFFT`, BART's `-w` option of the commands, converted to the tolerance that gives that width | -| Oversampling $\sigma$ | 1.25; 2 for `ncalib`, `nlinv` and `rtnlinv` | `oversampling` of {class}`~bartorch.linop.NUFFT`, BART's `-o` option | +| Tolerance $\varepsilon$ | $10^{-3}$; $10^{-6}$ for `ncalib` and `nlinv` | The kernel width: `width` of {class}`~bartorch.linop.NUFFT`, BART's `-w` option of the commands, converted to the tolerance that gives that width | +| Oversampling $\sigma$ | 1.25; 2 for `ncalib` and `nlinv` | `oversampling` of {class}`~bartorch.linop.NUFFT`, BART's `-o` option | The relative error of a transform is of the order of $\varepsilon$; the test suite compares it with the explicit sum above. The default trades accuracy diff --git a/docs/guides/user/installation.md b/docs/guides/user/installation.md index 086afcee..0511d429 100644 --- a/docs/guides/user/installation.md +++ b/docs/guides/user/installation.md @@ -5,7 +5,7 @@ Install the PyTorch build for the target device first, with the command the [PyTorch installation selector](https://pytorch.org/get-started/locally/) gives for the CPU or for a CUDA version. bartorch requires PyTorch 2.2 or -later; installing bartorch into an environment without PyTorch installs the +later (2.3 on macOS, 2.7.1 on Linux); installing bartorch into an environment without PyTorch installs the default build from PyPI. ## bartorch installation @@ -146,12 +146,13 @@ and the configurations that are refused. ## OpenMP runtime BART and FINUFFT are compiled with OpenMP and share one runtime with PyTorch. -On macOS and Windows the library is linked against the runtime PyTorch -contains (`torch/lib/libomp.dylib`, `torch/lib/libiomp5md.dll`) and contains -none of its own: the LLVM runtime terminates the process when a second copy +The library is linked against the runtime PyTorch contains +(`torch/lib/libgomp.so.1`, `torch/lib/libomp.dylib`, +`torch/lib/libiomp5md.dll`) and the wheels contain none of their own; on macOS +and Windows the LLVM runtime terminates the process when a second copy initializes (`OMP: Error #15`). No installed file is modified, and -`KMP_DUPLICATE_LIB_OK` is not used. On Linux the library uses the compiler's -runtime, `libgomp` for the GCC-built wheel. +`KMP_DUPLICATE_LIB_OK` is not used. On Linux this needs PyTorch 2.7.1 or +later, the first release whose copy is named `libgomp.so.1`. ## Command-line interface diff --git a/docs/guides/user/prerequisites.md b/docs/guides/user/prerequisites.md index cf6bc714..5555c636 100644 --- a/docs/guides/user/prerequisites.md +++ b/docs/guides/user/prerequisites.md @@ -5,7 +5,7 @@ | Requirement | Version | | --- | --- | | Python | 3.10 or later; each wheel is tested on 3.10 and 3.14 | -| PyTorch | 2.2 or later (2.3 on macOS, the first to carry `libomp.dylib`), CPU or CUDA build | +| PyTorch | 2.2 or later (2.3 on macOS, the first to carry `libomp.dylib`; 2.7.1 on Linux, the first whose `libgomp` is named `libgomp.so.1`), CPU or CUDA build | | NumPy, SciPy | NumPy 1.24 and SciPy 1.10 or later; MRI-NUFFT, a dependency, raises these to NumPy 2.2 and SciPy 1.13 | | TorchSim, MRI-NUFFT | TorchSim 0.0.5 and MRI-NUFFT 1.0 or later, installed as dependencies | diff --git a/docs/misc/license.md b/docs/misc/license.md index 9bf7bea6..3ae55099 100644 --- a/docs/misc/license.md +++ b/docs/misc/license.md @@ -16,7 +16,6 @@ replace; their notices are distributed with the source and in each wheel's | BART | `external/bart/` (Git submodule, `pulserver/bart`) | BSD-3-Clause, `external/bart/LICENSE`, with further notices in individual files | | pocketfft | `external/pocketfft/` | BSD-3-Clause, `external/pocketfft/LICENSE.md` | | BlocksRuntime (LLVM compiler-rt) | `external/blocksruntime/` | University of Illinois/NCSA or MIT, at the user's choice, `external/blocksruntime/LICENSE.TXT`; linked into Linux builds made with clang | -| GNU OpenMP runtime (`libgomp`) | Added to the Linux wheel by `auditwheel` | GPL-3.0 with the GCC Runtime Library Exception | BART's license: diff --git a/external/bart b/external/bart index 694ab5fe..7fab96d7 160000 --- a/external/bart +++ b/external/bart @@ -1 +1 @@ -Subproject commit 694ab5fe8241c43d67df5382c38dcf074aaa95a3 +Subproject commit 7fab96d7369e0be35c2146533476f59c3ba14d53 diff --git a/pyproject.toml b/pyproject.toml index 0fe3949c..5fd2d923 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,9 @@ dependencies = [ # On macOS the library's OpenMP is the libomp.dylib torch carries, which # torch's wheels carry from 2.3 on. "torch>=2.3; sys_platform == 'darwin'", + # On Linux it is torch's libgomp.so.1; before 2.7.1 torch's copy carried + # a hashed name the library's NEEDED entry does not match. + "torch>=2.7.1; sys_platform == 'linux'", # TorchSim is the signal model. `nlop.FromTorchSim` turns any of its # simulators into a BART nonlinear operator -- value, Jacobian-vector # product and adjoint product, which is exactly what `nlop_s` asks for -- diff --git a/scripts/gen_catalogue.py b/scripts/gen_catalogue.py index 330b6072..fbdbff01 100644 --- a/scripts/gen_catalogue.py +++ b/scripts/gen_catalogue.py @@ -60,7 +60,7 @@ class Unknown(Exception): TYPED_KINDS = ( "INT UINT LONG ULONG ULLONG PINT FLOAT DOUBLE STRING " "INFILE OUTFILE INOUTFILE CFL " - "VEC2 VEC3 VECN FLVEC2 FLVEC3 FLVEC4 FLVECN DOVEC3 DOVECN" + "VEC2 VEC3 VECN FLVEC2 FLVEC3 FLVEC4 FLVEC7 FLVECN DOVEC3 DOVECN" ).split() #: The two that take a variable number of values and so have no metavar of @@ -263,7 +263,7 @@ def read_command(name: str, source: Path) -> dict: """One BART tool, as its own source declares it.""" code = strip_comments(source.read_text(errors="replace")) - helped = re.search(r"help_str\[\]\s*=\s*((?:\s*\"(?:[^\"\\]|\\.)*\")+)", code) + helped = re.search(r"help_str(?:\[\])?\s*=\s*((?:\s*\"(?:[^\"\\]|\\.)*\")+)", code) arguments = [ read_argument(macro, args) for macro, args in macro_calls( diff --git a/src/bartorch/_call.py b/src/bartorch/_call.py index 3c00afe3..d1e3af2e 100644 --- a/src/bartorch/_call.py +++ b/src/bartorch/_call.py @@ -89,9 +89,8 @@ def convert(self, value: Any, inputs: list, command: str) -> Any: #: BART's term letters by parser: ``pics``, ``wshfl`` and ``denoise`` share -#: one; ``sqpics`` and ``moba`` each have their own. +#: one; ``moba`` has its own. _PICS_KINDS = frozenset("W H N L T G C V P S Q F I R1 R2".split()) -_SQPICS_KINDS = frozenset("W L T I Q R1 R2".split()) _MOBA_KINDS = frozenset("W T Q".split()) #: Each argument of a public command that BART takes as dimensions or a @@ -132,7 +131,6 @@ def convert(self, value: Any, inputs: list, command: str) -> Any: ("pics", "-L"): _Takes("axes", "Axes reconstructed one at a time (batch mode)."), ("pics", "--shared-img-dims"): _Takes("axes", "Axes the image is shared along."), ("pics", "--mpi"): _Takes("axes", "Axes distributed over MPI processes."), - ("sqpics", "-R"): _Takes("regularizers", "Regularization terms.", kinds=_SQPICS_KINDS), ("ssa", "-g"): _Takes("indices", "Grouping, as a set of indices."), ("wshfl", "-R"): _Takes("regularizers", "Regularization terms.", kinds=_PICS_KINDS), } @@ -185,6 +183,7 @@ def translate(command: str, values: dict[str, Any], inputs: list) -> dict[str, A "FLVEC2": "tuple[float, float]", "FLVEC3": "tuple[float, float, float]", "FLVEC4": "tuple[float, float, float, float]", + "FLVEC7": "tuple[float, float, float, float, float, float, float]", "FLVECN": "tuple[float, ...]", "DOVEC3": "tuple[float, float, float]", "DOVECN": "tuple[float, ...]", @@ -280,16 +279,17 @@ def signature_for(name: str) -> inspect.Signature: return inspect.Signature(parameters, return_annotation="torch.Tensor | tuple | str | None") -#: Commands whose only output BART marks optional and then creates anyway. +#: Commands whose only output BART marks optional and then needs anyway. # # `create_cfl` hands the name to `io_unlink_if_opened`, which calls `strcmp` # on it, so an omitted name is a segmentation fault rather than a tool that -# quietly writes nothing -- `mobafit.c:398` and `morphop.c:77` are both -# unguarded. Every other command with an optional output writes it through -# `anon_cfl` or behind an `if`, prints its answer instead when no name is -# given, and is left alone: passing one would silently turn the printed line -# `estdelay` and `measure` return here into a tensor. -_WRITES_ITS_OPTIONAL_OUTPUT = frozenset({"mobafit", "morphop"}) +# quietly writes nothing -- `morphop.c:77` is unguarded. `raga` refuses an +# omitted name (`raga.c:120`) unless `--search-tiny` asks it to print instead. +# Every other command with an optional output writes it through `anon_cfl` or +# behind an `if`, prints its answer instead when no name is given, and is left +# alone: passing one would silently turn the printed line `estdelay` and +# `measure` return here into a tensor. +_WRITES_ITS_OPTIONAL_OUTPUT = frozenset({"morphop", "raga"}) def _outputs(command: Command) -> int: @@ -394,7 +394,9 @@ def _docstring(command: Command, parameters: list[inspect.Parameter]) -> str: continue lines.append(f"{keyword} : {parameter.annotation}{_default(parameter)}") rule = TRANSLATED.get((command.name, option.flag)) - said = rule.help if rule is not None else option.help.strip() or f"BART's {option.flag}." + # `|M|` in BART's help would read as a reStructuredText substitution. + said = rule.help if rule is not None else option.help.strip().replace("|", r"\|") + said = said or f"BART's {option.flag}." lines.append(f" {said} (``{option.flag}``)") lines += ["**extra : Any", " Further BART flags, passed through by name."] diff --git a/src/bartorch/_catalogue.py b/src/bartorch/_catalogue.py index a16ac1aa..75948339 100644 --- a/src/bartorch/_catalogue.py +++ b/src/bartorch/_catalogue.py @@ -221,6 +221,60 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('T', '', 'SET', '', 'varying number of threads'), Option('S', '', 'SET', '', 'varying problem size'), Option('s', '', 'ULONG', 'flags', 'select benchmarks'), + Option('g', '', 'SET', '', 'perform benchmark on GPU'), + ), + ), + 'bet': Command( + name='bet', + help='Create 3D mask from 3D brain image using BET algorithm.', + arguments=( + Argument('brain image (complex or magnitude)', 'INFILE', True), + Argument('output', 'OUTFILE', True), + ), + options=( + Option( + 's', + '', + 'INT', + 'sub', + 'number of subdivisons for initialization of the polyedron - default 4 (needs to be in range of 2 - 5)', + ), + Option( + 'e', + '', + 'INT', + 'erode', + 'number or erosion steps - default 0 (needs to be in range of 0 - 10)', + ), + Option( + 'b', + '', + 'FLOAT', + 'bt', + 'fractional intensity threshold; between 0 and 1 - default 0.5', + ), + Option( + '', + 'res', + 'FLVEC3', + 'res_x:res_y:res_z', + 'resolution/voxel size [mm] - default 1.0:1.0:1.0', + ), + Option( + 'i', + '', + 'INT', + 'n_iter', + 'number of iterations for the BET algorithm - default 1000', + ), + Option('', 'stl_out', 'OUTFILE', 'stl_outfile', 'export stl mesh to file'), + Option( + 'w', + '', + 'SET', + '', + 'use winding number method for mesh to mask conversion (default: false)', + ), ), ), 'bin': Command( @@ -296,11 +350,19 @@ def choices(self) -> dict[str, tuple[Option, ...]]: help='Perform function evaluation on array.', arguments=( Argument('func', 'STRING', True), - Argument('input', 'INFILE', True), + Argument('input1', 'INFILE', True), + Argument('input2', 'INFILE', False), Argument('output', 'OUTFILE', True), ), options=( Option('L', '', 'SPECIAL', '', 'Print a list of all supported functions'), + Option( + 'v', + '', + 'CFL', + 'scalar', + 'scalar value for functions that require a scalar argument', + ), ), ), 'caldir': Command( @@ -384,6 +446,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: 'dim', 'Perform alignment of coil sensitivities along dimension A', ), + Option('W', '', 'INFILE', '', 'pre-whiten k-space before cc'), ), ), 'cdf97': Command( @@ -841,6 +904,15 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Argument('output', 'OUTFILE', True), ), ), + 'extractdc': Command( + name='extractdc', + help='Extract DC component from radial k-space data correcting for shifted or asymmetric trajectories.', + arguments=( + Argument('traj', 'INFILE', True), + Argument('kspace', 'INFILE', True), + Argument('dc', 'OUTFILE', True), + ), + ), 'fakeksp': Command( name='fakeksp', help='Recreate k-space from image and sensitivities.', @@ -1165,6 +1237,24 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '&(config.slice_ord)', 'ISMRMRD_SLICE_INTERLEAVED_SIEMENS', ), + Option( + '\\0', + 'write-ismrm-cfl-img', + 'SELECT', + '', + 'Output ISMRM data as cfl image.', + '&write_ismrm', + 'ISMRM_WRITE_CFL_IMG', + ), + Option( + '\\0', + 'write-ismrm-ush-img', + 'SELECT', + '', + 'Output ISMRM data as unsigned short image.', + '&write_ismrm', + 'ISMRM_WRITE_USHORT_IMG', + ), Option('\\0', 'stream', 'SET', '', 'Use streaming protocols.'), ), ), @@ -1412,7 +1502,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '', 'UINT', 'model', - 'Select the MGRE model from enum { WF = 0, WFR2S, WF2R2S, R2S, PHASEDIFF, ..., WF_fB0, WF_R2S, T1_R2S, W_T1_F_T1_RS2 } [default: WFR2S]', + 'Select the MGRE model from enum { WF = 0, WFR2S, WF2R2S, R2S, PHASEDIFF, PI, T1_R2S, W_T1_F_T1_RS2 } [default: WFR2S]', ), Option('l', '', 'PINT', '1/-l2', 'toggle l1-wavelet or l2 regularization.'), Option('i', '', 'PINT', 'iter', 'Number of Newton steps'), @@ -1420,21 +1510,29 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('j', '', 'FLOAT', 'minreg', 'Minimum regularization parameter'), Option('u', '', 'FLOAT', 'rho', 'ADMM rho [default: 0.01]'), Option('C', '', 'PINT', 'iter', 'inner iterations'), + Option('', 'tol', 'FLOAT', 'tol', 'tolerance for fista early stopping (default: 0.01)'), Option('s', '', 'FLOAT', 'step', 'step size'), Option('B', '', 'FLOAT', 'bound', 'lower bound for relaxation'), Option( 'b', '', 'FLVEC2', - 'SMO:SC', - 'B0 field: spatial smooth level; scaling [default: 222.; 1.]', + 'a:b', + 'B0 field: sobolev parameter (a=0 means no sobolev) [default: 222.; 32.]', ), Option('d', '', 'INT', 'level', 'Debug level'), Option('N', '', 'SET', '', '(normalize)'), Option('f', '', 'FLOAT', 'FOV', ''), Option('p', '', 'INFILE', 'PSF', ''), Option('J', '', 'SET', '', 'Stack frames for joint recon'), - Option('M', '', 'SET', '', 'Simultaneous Multi-Slice reconstruction'), + Option('M', 'sms', 'SET', '', 'Simultaneous Multi-Slice reconstruction'), + Option( + '', + 'sos', + 'SET', + '', + 'SoS reconstruction (FFTuc and Sobolev in slice dimension)', + ), Option( 'O', '', @@ -1443,7 +1541,13 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '(Output original maps from reconstruction without post processing)', ), Option('g', '', 'SET', '', 'use gpu'), - Option('', 'positive-maps', 'LONG', 'flag', 'Maps with positivity constraint as FLAG!'), + Option( + '', + 'positive-maps', + 'ULONG', + 'flag', + 'Maps with positivity constraint as FLAG!', + ), Option( '', 'not-wav-maps', @@ -1465,8 +1569,9 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('', 'multi-gpu', 'INT', 'num', '(number of gpus to use)'), Option('I', '', 'INFILE', 'init', 'File for initialization'), Option('t', '', 'INFILE', 'traj', 'K-space trajectory'), + Option('', 'sens', 'INFILE', 'sens', 'Use precomputed sensitivities'), Option('o', '', 'FLOAT', 'os', 'Oversampling factor for gridding [default: 1.]'), - Option('', 'img_dims', 'VEC3', 'x:y:z', 'dimensions'), + Option('x', 'img_dims', 'VEC3', 'x:y:z', 'dimensions'), Option('k', '', 'SET', '', 'k-space edge filter for non-Cartesian trajectories'), Option( '', @@ -1518,12 +1623,20 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '&conf.fat_spec', 'FAT_SPEC_0', ), + Option( + '', + 'field-strength', + 'FLOAT', + 'B0', + 'B0 field strength for fat spectrum [def: 3T]', + ), Option('', 'scale_data', 'FLOAT', '', 'scaling factor for data'), Option('', 'scale_psf', 'FLOAT', '', '(scaling factor for PSF)'), Option('', 'normalize_scaling', 'SET', '', '(normalize scaling by data / PSF)'), Option('', 'seq', 'SUBOPT', '...', 'configure sequence parameters'), Option('', 'sim', 'SUBOPT', '...', 'configure simulation parameters'), Option('', 'other', 'SUBOPT', '...', 'configure other parameters'), + Option('', 'monitor', 'OUTFILE', 'monitor', 'File for monitoring output'), ), ), 'mobafit': Command( @@ -1532,7 +1645,8 @@ def choices(self) -> dict[str, tuple[Option, ...]]: arguments=( Argument('enc', 'INFILE', True), Argument('echo/contrast images', 'INFILE', True), - Argument('coefficients', 'OUTFILE', False), + Argument('coefficients', 'OUTFILE', True), + Argument('covariance matrix', 'OUTFILE', False), ), options=( Option('B', '', 'SELECT', '', 'bSSFP', '&seq', 'BSSFP'), @@ -1558,6 +1672,15 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('G', '', 'SELECT', '', 'MGRE', '&seq', 'MGRE'), Option('D', '', 'SELECT', '', 'diffusion', '&seq', 'DIFF'), Option('S', '', 'SELECT', '', 'Simulation based fitting', '&seq', 'SIM'), + Option( + 'P', + '', + 'SELECT', + '', + 'Phase fitting f(p0, fB0) = |M| exp(i 2pi p0) * exp(i 2pi * fB0 * te)', + '&seq', + 'PHASE', + ), Option( 'm', '', @@ -1583,6 +1706,13 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '', 'Initial values of parameters in model-based reconstruction', ), + Option( + '', + 'init-file', + 'INFILE', + 'init', + 'Initial values of parameters in model-based reconstruction', + ), Option('', 'scale', 'FLVECN', '', 'Scaling'), Option( '', @@ -1618,6 +1748,13 @@ def choices(self) -> dict[str, tuple[Option, ...]]: ), Option('', 'b1map', 'INFILE', '[deg]', 'Input B1 map as cfl file'), Option('', 'b0map', 'INFILE', '[rad/s]', 'Input B0 map as cfl file'), + Option( + '', + 'field-strength', + 'FLOAT', + 'B0', + 'B0 field strength for fat spectrum [def: 3T]', + ), Option( '', 'seq', @@ -1629,6 +1766,22 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('', 'other', 'SUBOPT', '...', 'configure other simulation parameters'), Option('', 'pool', 'SUBOPT', '...', 'configure pool parameters for BMC simulation'), Option('', 'cest', 'SUBOPT', '...', 'configure parameters for CEST simulation'), + Option( + '', + 'fat_spec_0', + 'SELECT', + '', + 'select fat spectrum from ISMRM fat-water tool', + '&fat_spec', + 'FAT_SPEC_0', + ), + Option( + '', + 'fB0-init', + 'SET', + '', + 'Take fB0 average of echos with equal dTE for initialization of fB0', + ), ), ), 'mobasig': Command( @@ -1677,6 +1830,22 @@ def choices(self) -> dict[str, tuple[Option, ...]]: 'model', 'Select the MGRE model from enum { WF = 0, WFR2S, WF2R2S, R2S, PHASEDIFF } [default: WFR2S]', ), + Option( + '', + 'field-strength', + 'FLOAT', + 'B0', + 'B0 field strength for fat spectrum [def: 3T]', + ), + Option( + '', + 'fat_spec_0', + 'SELECT', + '', + 'select fat spectrum from ISMRM fat-water tool', + '&(data.fat_spec)', + 'FAT_SPEC_0', + ), ), ), 'morphop': Command( @@ -1692,6 +1861,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('d', '', 'SELECT', '', 'DILATION', '&morph_type', 'DILATION'), Option('o', '', 'SELECT', '', 'OPENING', '&morph_type', 'OPENING'), Option('c', '', 'SELECT', '', 'CLOSING', '&morph_type', 'CLOSING'), + Option('s', '', 'SELECT', '', 'SKELETON', '&morph_type', 'SKELETON'), Option('l', '', 'SELECT', '', 'LABEL', '&morph_type', 'LABEL'), Option('B', '', 'SELECT', '', 'use BALL structuring element', '&mask_type', 'BALL'), ), @@ -1730,6 +1900,14 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('c', '', 'FLOAT', '', '(c in 1 + a * \\Laplace^-b/2)'), Option('w', '', 'FLOAT', '', '(inverse scaling of the data)'), Option('o', '', 'SET', '', 'return oversampled coils'), + Option('', 'sms', 'SET', '', 'sms reconstruction (FFT in slice dimension)'), + Option( + '', + 'sos', + 'SET', + '', + 'sos reconstruction (ucFFT and Sobolev regularization in slice dimension)', + ), Option('N', '', 'SET', '', 'Normalize coil sensitivities'), Option('m', '', 'PINT', 'nmaps', 'Number of ENLIVE maps to use in reconstruction'), Option('x', 'dims', 'VEC3', 'x:y:z', 'Explicitly specify sens dimensions'), @@ -1750,6 +1928,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: 'd', 'Use phase pole detection after d iterations (0 for every iteration)', ), + Option('', 'nufft-conf', 'SUBOPT', '...', 'configure nufft'), ), ), 'nlinv': Command( @@ -2298,23 +2477,8 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('', 'ist_last', 'SET', '', 'end iteration with call to data consistency'), Option('', 'motion-field', 'INFILE', 'file', 'motion field'), Option('', 'nufft-conf', 'SUBOPT', '...', 'configure nufft'), - ), - ), - 'pocsense': Command( - name='pocsense', - help='Perform POCSENSE reconstruction.', - arguments=( - Argument('kspace', 'INFILE', True), - Argument('sensitivities', 'INFILE', True), - Argument('output', 'OUTFILE', True), - ), - options=( - Option('i', '', 'INT', 'iter', 'max. number of iterations'), - Option('r', '', 'FLOAT', 'alpha', 'regularization parameter'), - Option('l', '', 'PINT', '1/-l2', 'toggle l1-wavelet or l2 regularization'), - Option('g', '', 'SET', '', '()'), - Option('o', '', 'FLOAT', '', '()'), - Option('m', '', 'FLOAT', '', '()'), + Option('', 'field-map', 'INFILE', 'file', 'field map'), + Option('', 'time-map', 'INFILE', 'file', 'time map'), ), ), 'poisson': Command( @@ -2415,12 +2579,27 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('', 'N', 'PINT', 'int', 'number of time-steps (default = 1e6 * dur)'), ), ), + 'pulseq': Command( + name='pulseq', + help='Read/write pulseq files.', + arguments=( + Argument('pulseq-file', 'INOUTFILE', True), + Argument('gradients', 'OUTFILE', False), + Argument('adcs', 'OUTFILE', False), + Argument('pulse', 'OUTFILE', False), + ), + options=( + Option('w', '', 'SET', '', 'write'), + Option('F', '', 'INOUTFILE', 'file', 'RF Shapes file'), + Option('E', '', 'INOUTFILE', 'file', 'Events file'), + ), + ), 'raga': Command( name='raga', help='Generate file with RAGA indices for given approximated tiny golden ratio angle/raga increment and full frame spokes.', arguments=( Argument('spokes', 'INT', True), - Argument('output', 'OUTFILE', True), + Argument('output', 'OUTFILE', False), ), options=( Option('s', 'tiny-angle', 'PINT', '# Tiny GA', 'tiny (small) golden ratio angle'), @@ -2431,6 +2610,13 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('i', 'inversions', 'LONG', 'i', 'Number of (non-aligned) inversions'), Option('c', 'shifts', 'LONG', 'c', 'Number of (non-aligned) shifts'), Option('', 'no-double-base', 'CLEAR', '', 'Define GA over Pi base instead of 2Pi.'), + Option( + '', + 'search-tiny', + 'SET', + '', + 'Search for tiny golden-angle index with given number of spokes.', + ), ), ), 'reconet': Command( @@ -2540,7 +2726,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: ), 'reshape': Command( name='reshape', - help='Reshape selected dimensions.', + help='Reshape selected dimensions. One dimension can be set to -1 or 0, which will be filled in to match the total size.', arguments=( Argument('flags', 'ULONG', True), Argument('dim', 'TUPLE', True), @@ -2623,45 +2809,6 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Argument('output', 'OUTFILE', True), ), ), - 'rtnlinv': Command( - name='rtnlinv', - help='Jointly estimate a time-series of images and sensitivities with nonlinear\ninversion using {iter} iteration steps. Optionally outputs\nthe sensitivities.', - arguments=( - Argument('kspace', 'INFILE', True), - Argument('output', 'OUTFILE', True), - Argument('sensitivities', 'OUTFILE', False), - ), - options=( - Option('i', '', 'UINT', 'iter', 'Number of Newton steps'), - Option('R', '', 'FLOAT', '', '(reduction factor)'), - Option('M', '', 'FLOAT', '', '(minimum for regularization)'), - Option('d', '', 'INT', 'level', 'Debug level'), - Option('c', '', 'SET', '', 'Real-value constraint'), - Option('N', '', 'CLEAR', '', 'Do not normalize image with coil sensitivities'), - Option('m', '', 'PINT', 'nmaps', 'Number of ENLIVE maps to use in reconstruction'), - Option('U', '', 'CLEAR', '', 'Do not combine ENLIVE maps in output'), - Option('f', '', 'FLOAT', 'FOV', 'restrict FOV'), - Option('p', '', 'INFILE', 'file', 'pattern / transfer function'), - Option('t', '', 'INFILE', 'file', 'kspace trajectory'), - Option('I', '', 'INFILE', 'file', 'File for initialization'), - Option( - 'C', - '', - 'INFILE', - '', - '(File for initialization with image space sensitivities)', - ), - Option('g', '', 'SET', '', 'use gpu'), - Option('S', '', 'SET', '', 'Re-scale image after reconstruction'), - Option('a', '', 'FLOAT', '', '(a in 1 + a * \\Laplace^-b/2)'), - Option('b', '', 'FLOAT', '', '(b in 1 + a * \\Laplace^-b/2)'), - Option('T', '', 'FLOAT', 'temp_damp', 'temporal damping [default: 0.9]'), - Option('w', '', 'FLOAT', '', '(inverse scaling of the data)'), - Option('x', '', 'VEC3', 'x:y:z', 'Explicitly specify image dimensions'), - Option('A', '', 'SET', '', '(Alternative scaling)'), - Option('s', '', 'SET', '', '(Simultaneous Multi-Slice reconstruction)'), - ), - ), 'sake': Command( name='sake', help='Use SAKE algorithm to recover a full k-space from undersampled\ndata using low-rank matrix completion.', @@ -2738,8 +2885,9 @@ def choices(self) -> dict[str, tuple[Option, ...]]: ), 'seq': Command( name='seq', - help='Computes a GRE sequence.', + help='Computes a Sequence.', arguments=( + Argument('seq', 'STRING', False), Argument( '0th moment (x,y,z) at sample points, sample_points, phase of adc', 'OUTFILE', @@ -2749,176 +2897,6 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Argument('0th moment (x,y,z) per imaging block', 'OUTFILE', False), Argument('pulseq file', 'OUTFILE', False), ), - options=( - Option( - 'd', - '', - 'FLOAT', - 'dt', - 'time-increment per sample (default: seq->conf->phys.tr / 1000)', - ), - Option('N', '', 'LONG', 'samples', 'Number of samples (default: 1000)'), - Option('s', '', 'DOVEC3', 'RO:PE:SL', 'FOV shift'), - Option('S', '', 'DOVEC3', 'RO:PE:SL', 'relative FOV shift'), - Option( - '', - 'dist', - 'FLOAT', - 'dist', - 'slice distance factor [1 / slice_thickness] (default: 1.)', - ), - Option( - '', - 'no-spoiling', - 'SELECT', - '', - 'spoiling off (default: rf random)', - '&seq->conf->phys.contrast', - 'SEQ_CONTRAST_NO_SPOILING', - ), - Option( - '', - 'spoiled', - 'SELECT', - '', - 'RF_SPOILED (inc: 50 deg, gradient on) (default: rf random)', - '&seq->conf->phys.contrast', - 'SEQ_CONTRAST_RF_SPOILED', - ), - Option('', 'FOV', 'DOUBLE', 'FOV', 'Field Of View'), - Option('', 'BR', 'PINT', 'BR', 'Base Resolution'), - Option('', 'slice_thickness', 'DOUBLE', 'slice_thickness', 'Slice thickness'), - Option('', 'FA', 'DOUBLE', 'flip angle', 'Flip angle [deg]'), - Option('', 'TR', 'DOUBLE', 'TR', 'TR'), - Option('', 'TE', 'DOUBLE', 'TE', 'TE'), - Option('', 'TE_delta', 'DOUBLE', 'TE_delta', 'TE_delta'), - Option('', 'BWTP', 'DOUBLE', 'BWTP', 'Bandwidth Time Product'), - Option('', 'rf_duration', 'DOUBLE', 'rf_duration', 'RF pulse duration'), - Option('', 'dwell', 'DOUBLE', 'dwell', 'Dwell time'), - Option('', 'os', 'DOUBLE', 'os', 'Oversampling factor'), - Option('', 'pe_mode', 'UINT', 'pe_mode', 'Phase-encoding mode'), - Option( - '', - 'turn', - 'SELECT', - '', - 'turn-based PE (default: RAGA)', - '&seq->conf->enc.pe_mode', - 'SEQ_PEMODE_TURN', - ), - Option( - '', - 'mems', - 'SELECT', - '', - 'multi-echo/multi-spoke PE (default: RAGA)', - '&seq->conf->enc.pe_mode', - 'SEQ_PEMODE_MEMS_HYB', - ), - Option( - '', - 'raga', - 'SELECT', - '', - 'RAGA PE', - '&seq->conf->enc.pe_mode', - 'SEQ_PEMODE_RAGA', - ), - Option( - '', - 'raga_flags', - 'ULONG', - 'raga_aligned_flags', - 'RAGA aligned flags (by bitmask)', - ), - Option( - '', - 'chrono', - 'SET', - '', - 'save gradients/moments/sampling in chronological order (RAGA)', - ), - Option('R', '', 'OUTFILE', 'file', 'raga indices'), - Option('', 'tiny', 'PINT', 'tiny', 'Tiny golden-ratio index'), - Option('e', 'echoes', 'LONG', 'echoes', 'Number of echoes'), - Option('r', 'lines', 'LONG', 'lines', 'Number of phase encoding lines'), - Option( - 'z', - 'partitions', - 'LONG', - 'partitions', - 'Number of partitions (3D) or SMS groups (2D)', - ), - Option( - 't', - 'measurements', - 'LONG', - 'measurements', - 'Number of measurements / frames (RAGA: total number of spokes)', - ), - Option( - 'f', - 'raga_full_frames', - 'LONG', - 'raga_full_frames', - 'Number of full frames (only RAGA)', - ), - Option('m', 'slices', 'LONG', 'slices', 'Number of slices of multiband factor (SMS)'), - Option('i', 'inversions', 'LONG', 'inversions', 'Number of inversions'), - Option( - '', - 'sequential-multislice', - 'SELECT', - '', - 'seq_order: sequential multislice (default: avg outer)', - '&seq->conf->enc.order', - 'SEQ_ORDER_SEQ_MS', - ), - Option( - '', - 'avg-inner', - 'SELECT', - '', - 'seq_order: average inner (default: avg outer)', - '&seq->conf->enc.order', - 'SEQ_ORDER_AVG_INNER', - ), - Option('', 'mb_factor', 'PINT', 'mb_factor', 'Multi-band factor'), - Option('', 'sms_distance', 'DOUBLE', 'sms_distance', 'SMS slice distance'), - Option( - '', - 'IR_NON', - 'SELECT', - '', - 'Magn. preparation: Nonselective Inversion (default: off)', - '&seq->conf->magn.mag_prep', - 'SEQ_PREP_IR_NONSELECTIVE', - ), - Option('', 'TI', 'DOUBLE', 'TI', 'Inversion time'), - Option('', 'init_delay', 'DOUBLE', 'init_delay', 'Initial delay of measurement'), - Option('', 'inv_delay', 'DOUBLE', 'inv_delay_time', 'Inversion delay time'), - Option( - '', - 'gradient-normal', - 'SELECT', - '', - 'Gradient normal mode (default: fast)', - '&gradient_mode', - 'GRAD_NORMAL', - ), - Option( - '', - 'gradient-whisper', - 'SELECT', - '', - 'Gradient whispher mode (default: fast)', - '&gradient_mode', - 'GRAD_WHISPER', - ), - Option('', 'CUSTOM_LONG', 'VECN', '', 'custom long parameters'), - Option('', 'CUSTOM_DOUBLE', 'DOVECN', '', 'custom double parameters'), - Option('', 'LOOP', 'VECN', '', 'sequence loop dimensions'), - ), ), 'show': Command( name='show', @@ -3073,54 +3051,22 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Argument('output', 'OUTFILE', True), ), ), - 'spow': Command( - name='spow', - help='Raise array to the power of {exponent}. The exponent can be a complex number.', + 'sort': Command( + name='sort', + help='Sorts an array by real part and returns permutation to be used with bin.', arguments=( - Argument('exponent', 'CFL', True), Argument('input', 'INFILE', True), Argument('output', 'OUTFILE', True), ), ), - 'sqpics': Command( - name='sqpics', - help='Parallel-imaging compressed-sensing reconstruction.', + 'spow': Command( + name='spow', + help='Raise array to the power of {exponent}. The exponent can be a complex number.', arguments=( - Argument('kspace', 'INFILE', True), - Argument('sensitivities', 'INFILE', True), + Argument('exponent', 'CFL', True), + Argument('input', 'INFILE', True), Argument('output', 'OUTFILE', True), ), - options=( - Option('l', '', 'SPECIAL', '1/-l2', 'toggle l1-wavelet or l2 regularization.'), - Option('r', '', 'FLOAT', 'lambda', 'regularization parameter'), - Option( - 'R', - '', - 'SPECIAL', - ':A:B:C', - 'generalized regularization options (-Rh for help)', - ), - Option('s', '', 'FLOAT', 'step', 'iteration stepsize'), - Option('i', '', 'PINT', 'iter', 'max. number of iterations'), - Option('t', '', 'INFILE', 'file', 'k-space trajectory'), - Option('n', '', 'CLEAR', '', 'disable random wavelet cycle spinning'), - Option('g', '', 'SET', '', 'use GPU'), - Option('p', '', 'INFILE', 'file', 'pattern or weights'), - Option('I', '', 'SELECT', '', '(select IST)', '&ropts.algo', 'IST'), - Option('b', '', 'PINT', 'blk', 'Lowrank block size'), - Option('e', '', 'SET', '', 'Scale stepsize based on max. eigenvalue'), - Option('H', '', 'SET', '', '(hogwild)'), - Option('F', '', 'SET', '', '(fast)'), - Option('T', '', 'INFILE', 'file', '(truth file)'), - Option('W', '', 'INFILE', '', 'Warm start with '), - Option('d', '', 'INT', 'level', 'Debug level'), - Option('u', '', 'FLOAT', 'rho', 'ADMM rho'), - Option('C', '', 'PINT', 'iter', 'ADMM max. CG iterations'), - Option('f', '', 'FLOAT', 'rfov', 'restrict FOV'), - Option('m', '', 'SELECT', '', 'Select ADMM', '&ropts.algo', 'ADMM'), - Option('w', '', 'FLOAT', 'val', 'scaling'), - Option('S', '', 'SET', '', 'Re-scale the image after reconstruction'), - ), ), 'squeeze': Command( name='squeeze', @@ -3167,7 +3113,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: name='stl', help="Read and write stl files with '.stl' or cfl fileformat.", arguments=( - Argument('output', 'OUTFILE', True), + Argument('output', 'OUTFILE', False), ), options=( Option('', 'input', 'INFILE', '', 'Path to input file (.stl or cfl file format).'), @@ -3178,6 +3124,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: '', 'Internal stl model (help: bart stl --model h).', ), + Option('', 'cfov', 'SET', '', 'Scale and move model to center of FOV.'), Option( 's', '', @@ -3186,6 +3133,13 @@ def choices(self) -> dict[str, tuple[Option, ...]]: 'Multiplicate all coordinates of model with a scale factor.', ), Option('m', '', 'FLVEC3', 'move', 'Move model by vector.'), + Option( + '', + 'transform', + 'FLVEC7', + 'transform', + 'scale:move_x:move_y:move_z:rot_xy[deg]:rot_xz[deg]:rot_yz[deg]. Rotates the *centered* model.', + ), Option('', 'stat', 'SET', '', 'Show statistics of model.'), Option('', 'binary', 'CLEAR', '', 'Output STL files in binary format.'), Option( @@ -3436,6 +3390,7 @@ def choices(self) -> dict[str, tuple[Option, ...]]: Option('L', '', 'SET', '', 'use linectr offset'), Option('P', '', 'SET', '', 'use partctr offset'), Option('N', '', 'SET', '', 'only get noise'), + Option('D', '', 'SET', '', 'only get dummy (preparation scans)'), Option('R', '', 'SET', '', 'get data of reference scan'), Option('S', '', 'CLEAR', '', "don't include reference lines"), Option( @@ -3446,23 +3401,23 @@ def choices(self) -> dict[str, tuple[Option, ...]]: 'ignore (squash) selected dimensions (defaults to LEVEL_FLAG)', ), Option('C', '', 'SET', '', 'read data chronologically and ignore adc postitions'), - Option('', 'rational', 'SET', '', 'Rational Approximation Sampling'), Option('M', '', 'SET', '', 'MPI mode'), - Option('', 'bin', 'PINT', 'd', 'Binning of spokes for RAGA sampled data'), Option('X', '', 'CLEAR', '', 'no consistency check for number of read acquisitions'), Option('d', '', 'INT', 'level', 'Debug level'), + Option('', 'pattern', 'OUTFILE', 'pattern file', 'pattern file'), ), ), 'unwrap': Command( name='unwrap', help='Unwrap along selected dimensions.', arguments=( - Argument('dim', 'INT', True), + Argument('flags', 'ULONG', True), Argument('input', 'INFILE', True), Argument('output', 'OUTFILE', True), ), options=( Option('b', '', 'FLOAT', 'bounds', 'bounds (default: PI)'), + Option('l', '', 'SELECT', '', 'select Laplacian-based unwrapping', '&mode', 'MODE_LAP'), ), ), 'upat': Command( diff --git a/src/bartorch/_coverage.py b/src/bartorch/_coverage.py index bcd5dc19..2cab88d9 100644 --- a/src/bartorch/_coverage.py +++ b/src/bartorch/_coverage.py @@ -58,7 +58,9 @@ PRIVATE: dict[str, str] = { "bart": "the dispatcher that runs the other commands", "ismrmrd": "needs libismrmrd, which this build does not compile", - **dict.fromkeys(("tee", "multicfl", "tensorflow", "twixread", "toimg", "toraw"), _FILES), + **dict.fromkeys( + ("tee", "multicfl", "tensorflow", "twixread", "toimg", "toraw", "pulseq"), _FILES + ), **dict.fromkeys(("stl", "pol2mask", "morphop"), "mesh and mask geometry, not needed here"), **dict.fromkeys( ("cunet", "mnist", "nnet", "reconet", "nlinvnet", "sample", "onehotenc"), _NETWORK @@ -99,11 +101,12 @@ "svd", "hist", "compress", + "sort", ), _TORCH, ), **dict.fromkeys(("bloch", "epg", "mobasig", "pulse", "seq", "sim", "signal"), _SIMULATED), - **dict.fromkeys(("fftrot", "gmm"), _LATER), + **dict.fromkeys(("fftrot", "gmm", "bet", "extractdc"), _LATER), "ictv": "fails for every input in this BART (ictv.c:97 reshapes the wrong side)", "version": "bartorch.bart_version() reports it", "bitmask": "converts bitmasks, which the Python API does not use", diff --git a/src/bartorch/_dispatch.py b/src/bartorch/_dispatch.py index 8da550c9..2159dbb2 100644 --- a/src/bartorch/_dispatch.py +++ b/src/bartorch/_dispatch.py @@ -298,7 +298,7 @@ def build_argv( # that: `nlinv` and its relatives fit them from a low-resolution image that # everything after them is built on, and a grid twice over at FINUFFT's own # tolerance is a few megabytes there. So they get it. -_CALIBRATES = frozenset({"ncalib", "nlinv", "rtnlinv"}) +_CALIBRATES = frozenset({"ncalib", "nlinv"}) _CAREFUL_TOLERANCE = 1e-6 _CAREFUL_UPSAMPLING = 2.0 diff --git a/src/bartorch/_lib.py b/src/bartorch/_lib.py index 9cab353e..027d1d51 100644 --- a/src/bartorch/_lib.py +++ b/src/bartorch/_lib.py @@ -35,12 +35,11 @@ def _library_names() -> list[str]: def _load(path: Path) -> ctypes.CDLL: - if sys.platform not in ("win32", "darwin"): - return ctypes.CDLL(str(path)) - # The library imports the OpenMP runtime torch carries -- libiomp5md.dll - # on Windows, @rpath/libomp.dylib on macOS -- rather than one of its own. - # Importing torch first makes the copy torch loaded the one the loader - # binds to; torch's directory is where it is found otherwise. + # The library imports the OpenMP runtime torch carries -- libgomp.so.1 on + # Linux, @rpath/libomp.dylib on macOS, libiomp5md.dll on Windows -- rather + # than one of its own. Importing torch first makes the copy torch loaded + # the one the loader binds to; torch's directory is where it is found + # otherwise. import torch if sys.platform == "win32": diff --git a/src/bartorch/apps/pocsense.py b/src/bartorch/apps/pocsense.py index 5ec4e1c4..4cf58ddb 100644 --- a/src/bartorch/apps/pocsense.py +++ b/src/bartorch/apps/pocsense.py @@ -1,4 +1,4 @@ -"""``pocsense`` as a pipeline: the three projections, and the sweep over them.""" +"""POCSENSE as a pipeline: the three projections, and the sweep over them.""" from __future__ import annotations @@ -15,7 +15,7 @@ __all__ = ["pocsense"] -#: BART's own ``-i`` default (pocsense.c:66). +#: The sweep count BART's POCSENSE application defaults to. _MAXITER = 50 @@ -101,7 +101,7 @@ def pocsense( ) -> torch.Tensor: """POCSENSE reconstruction: the measured samples, the coils, and sparsity. - The pipeline :func:`bartorch.tools.pocsense` runs, assembled here: the + The pipeline of BART's POCSENSE application, assembled here: the scaling the application estimates, the sampling pattern read off the k-space, the modulation into the convention BART iterates in, and then a sweep of three projections from :mod:`bartorch.linop` and @@ -122,7 +122,7 @@ def pocsense( projections onto the data and onto the coils. wavelet : bool, default=False Threshold the wavelet coefficients of the coil images rather than - shrinking the samples towards zero, which is ``pocsense -l 1``. + shrinking the samples towards zero, which is the application's ``-l 1``. robust : float, default=None Soft-threshold the residual at a measured position by this much rather than discarding it, which is the application's ``-o``. diff --git a/src/bartorch/cli/_apps.py b/src/bartorch/cli/_apps.py index b3e38f6b..2fbfa22f 100644 --- a/src/bartorch/cli/_apps.py +++ b/src/bartorch/cli/_apps.py @@ -77,30 +77,8 @@ def _pics(options: dict[str, list[Any]], ndim: int) -> dict[str, Any]: return made -def _pocsense(options: dict[str, list[Any]], ndim: int) -> dict[str, Any]: - """``pocsense``'s flags as :func:`bartorch.apps.pocsense` takes them.""" - made: dict[str, Any] = {} - if "i" in options: - made["maxiter"] = _only(options.pop("i")) - if "r" in options: - made["alpha"] = _only(options.pop("r")) - if "o" in options: - made["robust"] = _only(options.pop("o")) - if "l" in options: - kind = _only(options.pop("l")) - if kind not in (1, 2): - raise Unsupported(f"pocsense -l takes 1 or 2, got {kind}") - made["wavelet"] = kind == 1 - if "g" in options: - options.pop("g") - if options: - raise Unsupported(f"pocsense {sorted(options)} is not something the app takes") - return made - - #: The commands an app answers, and the reader that fills its arguments. A #: command absent from here is BART's own. ADAPTERS = { "pics": _pics, - "pocsense": _pocsense, } diff --git a/src/bartorch/tools/recon.py b/src/bartorch/tools/recon.py index 0fb49a07..14f65d6e 100644 --- a/src/bartorch/tools/recon.py +++ b/src/bartorch/tools/recon.py @@ -258,10 +258,7 @@ def nlinv( "looklocker", "moba", "mobafit", - "pocsense", - "rtnlinv", "sake", - "sqpics", "wave", "wshfl", ) diff --git a/src/bartorch/util.py b/src/bartorch/util.py index 297c1d9a..cedcba35 100644 --- a/src/bartorch/util.py +++ b/src/bartorch/util.py @@ -308,7 +308,7 @@ def unwrap(input: torch.Tensor, axis: int, *, bound: float = math.pi) -> torch.T Half the period. """ dim = _bart_dim(axis, input.ndim) - return dispatch("unwrap", [input], None, _pos=[dim], b=float(bound)).reshape(input.shape) + return dispatch("unwrap", [input], None, _pos=[1 << dim], b=float(bound)).reshape(input.shape) @curated("casorati") diff --git a/src/csrc/abi/api.c b/src/csrc/abi/api.c index e01b7efb..8f5bab11 100644 --- a/src/csrc/abi/api.c +++ b/src/csrc/abi/api.c @@ -7,6 +7,8 @@ #include #include +#include + #include "misc/debug.h" #include "misc/misc.h" #include "misc/version.h" diff --git a/src/csrc/ops/ops.c b/src/csrc/ops/ops.c index dee19392..6539f73b 100644 --- a/src/csrc/ops/ops.c +++ b/src/csrc/ops/ops.c @@ -342,7 +342,7 @@ static int linop_nufft_worker(void* p) md_copy_dims(a->N, wgh_dims, a->wgh_dims); a->result = wrap_linop(nufft_create2(a->N, a->ksp_dims, a->cim_dims, a->traj_dims, a->traj, - wgh_dims, a->weights, a->bas_dims, a->basis, conf)); + wgh_dims, a->weights, a->bas_dims, a->basis, NULL, NULL, NULL, NULL, conf)); return 0; } diff --git a/src/csrc/ops/sense.c b/src/csrc/ops/sense.c index 58c16322..7e4edd3b 100644 --- a/src/csrc/ops/sense.c +++ b/src/csrc/ops/sense.c @@ -54,6 +54,7 @@ #include "noncart/nufft.h" #include "sense/model.h" +#include "sense/modelnc.h" #include "include/bartorch.h" @@ -95,6 +96,8 @@ extern const struct linop_s* bart_sense_nc_init(const bart_dim_t max_dims[DIMS], const bart_dim_t traj_dims[DIMS], const complex float* traj, const struct nufft_conf_s* conf, const bart_dim_t wgs_dims[DIMS], const complex float* weights, const bart_dim_t basis_dims[DIMS], const complex float* basis, + const bart_dim_t fieldmap_dims[DIMS], const complex float* fieldmap, + const bart_dim_t timemap_dims[DIMS], const complex float* timemap, const struct linop_s** fft_opp, bart_flags_t shared_img_dims); /* How many coils a slab holds. Zero leaves the operators as BART builds @@ -1440,6 +1443,8 @@ const struct linop_s* sense_nc_init(const bart_dim_t max_dims[DIMS], const bart_ const bart_dim_t traj_dims[DIMS], const complex float* traj, const struct nufft_conf_s* _conf, const bart_dim_t wgs_dims[DIMS], const complex float* weights, const bart_dim_t basis_dims[DIMS], const complex float* basis, + const bart_dim_t fieldmap_dims[DIMS], const complex float* fieldmap, + const bart_dim_t timemap_dims[DIMS], const complex float* timemap, const struct linop_s** fft_opp, bart_flags_t shared_img_dims) { bart_dim_t ksp_dims2[DIMS]; @@ -1448,13 +1453,15 @@ const struct linop_s* sense_nc_init(const bart_dim_t max_dims[DIMS], const bart_ bool sliced = sliceable((bart_dim_t)coil_batch, max_dims, map_dims, ksp_dims2, shared_img_dims) && ((NULL == weights) || (1 == wgs_dims[COIL_DIM])) - && ((NULL == basis) || (1 == basis_dims[COIL_DIM])); + && ((NULL == basis) || (1 == basis_dims[COIL_DIM])) + && (NULL == fieldmap) && (NULL == timemap) && !_conf->dft; if (!sliced) { chained(); return bart_sense_nc_init(max_dims, map_dims, maps, ksp_dims, traj_dims, traj, _conf, - wgs_dims, weights, basis_dims, basis, fft_opp, shared_img_dims); + wgs_dims, weights, basis_dims, basis, fieldmap_dims, fieldmap, + timemap_dims, timemap, fft_opp, shared_img_dims); } struct sense_s* d = sense_slabs((bart_dim_t)coil_batch, 0 != fold_maps, max_dims, map_dims, ksp_dims2, shared_img_dims, false); @@ -1466,7 +1473,7 @@ const struct linop_s* sense_nc_init(const bart_dim_t max_dims[DIMS], const bart_ d->maps = maps; d->slab = nufft_create2(DIMS, slab_ksp_dims, d->cim_dims, traj_dims, traj, (weights ? wgs_dims : NULL), weights, - (basis ? basis_dims : NULL), basis, *_conf); + (basis ? basis_dims : NULL), basis, NULL, NULL, NULL, NULL, *_conf); sense_output_from(d, false, 0, NULL); @@ -1675,7 +1682,7 @@ static const struct linop_s* form_transform(const struct bartorch_encoding* f, c return nufft_create2(DIMS, f->ksp_dims, cim_dims, f->traj_dims, f->traj, (f->weights ? f->wgh_dims : NULL), f->weights, - (f->basis ? f->bas_dims : NULL), f->basis, *conf); + (f->basis ? f->bas_dims : NULL), f->basis, NULL, NULL, NULL, NULL, *conf); } error("bartorch: %d is not one of this library's encoding transforms\n", f->transform); @@ -1698,7 +1705,7 @@ static const struct linop_s* form_chain(const struct bartorch_encoding* f, const if (BARTORCH_ENCODING_NUFFT == f->transform) return bart_sense_nc_init(max_dims, map_dims, f->sens, f->ksp_dims, f->traj_dims, f->traj, conf, - f->wgh_dims, f->weights, f->bas_dims, f->basis, NULL, 0); + f->wgh_dims, f->weights, f->bas_dims, f->basis, NULL, NULL, NULL, NULL, NULL, 0); if (0 != f->modulated) { diff --git a/src/csrc/substitute/nufft_finufft.c b/src/csrc/substitute/nufft_finufft.c index d97089ba..9dc499c1 100644 --- a/src/csrc/substitute/nufft_finufft.c +++ b/src/csrc/substitute/nufft_finufft.c @@ -90,7 +90,7 @@ extern void bartorch_paired_back(const struct bartorch_paired* p, int k, #include "include/bartorch.h" -extern struct linop_s* bart_nufft_create2(int N, const bart_dim_t ksp_dims[N], const bart_dim_t cim_dims[N], const bart_dim_t traj_dims[N], const complex float* traj, const bart_dim_t wgh_dims[N], const complex float* weights, const bart_dim_t bas_dims[N], const complex float* basis, struct nufft_conf_s conf); +extern struct linop_s* bart_nufft_create2(int N, const bart_dim_t ksp_dims[N], const bart_dim_t cim_dims[N], const bart_dim_t traj_dims[N], const complex float* traj, const bart_dim_t wgh_dims[N], const complex float* weights, const bart_dim_t bas_dims[N], const complex float* basis, const bart_dim_t fm_dims[N], const complex float* fieldmap, const bart_dim_t tm_dims[N], const complex float* timemap, struct nufft_conf_s conf); extern int bart_nufft_get_psf_dims(const struct linop_s* nufft, int N, bart_dim_t psf_dims[N]); extern void bart_nufft_get_psf(const struct linop_s* nufft, int N, const bart_dim_t psf_dims[N], complex float* psf); extern void bart_nufft_get_psf2(const struct linop_s* nufft, int N, const bart_dim_t psf_dims[N], const bart_stride_t psf_strs[N], complex float* psf); @@ -2595,7 +2595,7 @@ static const struct linop_s* toeplitz_for(int N, const bart_dim_t ksp_dims[N], c const struct linop_s* op = bart_nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, (NULL == basis) ? NULL : (along_samples ? stand_dims : bas_dims), - along_samples ? stand : basis, barts); + along_samples ? stand : basis, NULL, NULL, NULL, NULL, barts); struct nufft_data* data = CAST_DOWN(nufft_data, linop_get_data_nested(op)); @@ -2994,8 +2994,19 @@ static struct linop_s* try_create(int N, const bart_dim_t ksp_dims[N], const bar struct linop_s* nufft_create2(int N, const bart_dim_t ksp_dims[N], const bart_dim_t cim_dims[N], const bart_dim_t traj_dims[N], const complex float* traj, const bart_dim_t wgh_dims[N], const complex float* weights, - const bart_dim_t bas_dims[N], const complex float* basis, struct nufft_conf_s conf) + const bart_dim_t bas_dims[N], const complex float* basis, + const bart_dim_t fm_dims[N], const complex float* fieldmap, + const bart_dim_t tm_dims[N], const complex float* timemap, struct nufft_conf_s conf) { + /* `dft` is BART's explicit sum rather than its gridder, and the only + * transform that takes a field map. */ + if (conf.dft) + return bart_nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, + bas_dims, basis, fm_dims, fieldmap, tm_dims, timemap, conf); + + if ((NULL != fieldmap) || (NULL != timemap)) + error("bartorch: a field map needs the explicit transform (--nufft-conf dft).\n"); + struct linop_s* op = try_create(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, bas_dims, basis, conf); if (NULL != op) @@ -3011,7 +3022,7 @@ struct linop_s* nufft_create2(int N, const bart_dim_t ksp_dims[N], const bart_di count(CNT_BART); - return bart_nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, bas_dims, basis, barts_conf(conf)); + return bart_nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, bas_dims, basis, NULL, NULL, NULL, NULL, barts_conf(conf)); } struct linop_s* nufft_create(int N, const bart_dim_t ksp_dims[N], const bart_dim_t cim_dims[N], @@ -3021,7 +3032,7 @@ struct linop_s* nufft_create(int N, const bart_dim_t ksp_dims[N], const bart_dim bart_dim_t wgh_dims[N]; md_select_dims(N, ~MD_BIT(0), wgh_dims, traj_dims); - return nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, NULL, NULL, conf); + return nufft_create2(N, ksp_dims, cim_dims, traj_dims, traj, wgh_dims, weights, NULL, NULL, NULL, NULL, NULL, NULL, conf); } /* The rest read BART's own operator internals, so they are only safe on one diff --git a/src/csrc/substitute/psf.c b/src/csrc/substitute/psf.c index c65f006e..456d5b9f 100644 --- a/src/csrc/substitute/psf.c +++ b/src/csrc/substitute/psf.c @@ -170,7 +170,7 @@ static complex float* psf_int(int N, const bart_dim_t img_dims[N], const bart_di struct nufft_conf_s conf = psf_conf(periodic, lowmem, is_vptr(traj)); struct linop_s* op = nufft_create2(N, ksp_dims, img_dims2, trj_dims, traj, - wgh_dims, sqr_weights, sqr_bas_dims, sqr_basis, conf); + wgh_dims, sqr_weights, sqr_bas_dims, sqr_basis, NULL, NULL, NULL, NULL, conf); op = linop_reshape_in_F(op, N, img_dims); @@ -387,7 +387,7 @@ static complex float* psf_decomposed(bool to_host, bool real_out, const struct p * phase. Setting the points is the costly step: it sorts every * sample. */ struct linop_s* op = nufft_create2(ND, ksp_dims2, one_dims, trj_dims3, traj, - NULL, NULL, NULL, NULL, conf); + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, conf); bart_stride_t kstrs[ND]; bart_stride_t tstrs[ND]; diff --git a/tests/test_abi.py b/tests/test_abi.py index e7cf534a..fd1b8744 100644 --- a/tests/test_abi.py +++ b/tests/test_abi.py @@ -121,3 +121,26 @@ def test_build_info_reports_the_platforms_data_model(): info = dict(item.split("=", 1) for item in bartorch.build_info().split(",")) assert int(info["long"]) == ctypes.sizeof(ctypes.c_long) assert info["long"] == ("4" if sys.platform == "win32" else "8") + + +@pytest.mark.skipif(sys.platform == "win32", reason="a DLL exports only what it names") +def test_the_library_exports_no_c_symbol_but_the_bartorch_abi(): + """A C name BART leaves visible is one another library in the process can bind.""" + import shutil + import subprocess + + from bartorch._lib import library_path + + nm = shutil.which("nm") + if nm is None: + pytest.skip("no nm to read the dynamic symbol table with") + flags = ["-gU"] if sys.platform == "darwin" else ["-D", "--defined-only"] + command = [nm, *flags, str(library_path())] + listed = subprocess.run(command, capture_output=True, text=True, check=True) + names = [line.split()[-1] for line in listed.stdout.splitlines() if line.strip()] + if sys.platform == "darwin": + names = [name[1:] for name in names if name.startswith("_")] + # Names with a leading underscore or dot are the toolchain's: C++ mangling, + # the linker's section markers, OpenMP's named critical sections. + stray = sorted(n for n in names if not n.startswith(("_", ".", "bartorch_"))) + assert stray == [] diff --git a/tests/test_apps.py b/tests/test_apps.py index f5c951f4..b3a22b72 100644 --- a/tests/test_apps.py +++ b/tests/test_apps.py @@ -296,64 +296,50 @@ def test_the_fit_is_taken_from_where_it_is_started(): # --- pocsense -------------------------------------------------------------- -#: What the application is given, and what the app is given for the same run. -_POCS_CONFIGURATIONS = [ - ("plain", {}, {}), - ("l2", {"r": 0.1}, {"alpha": 0.1}), - ("wavelet", {"r": 0.01, "l": 1}, {"alpha": 0.01, "wavelet": True}), - ("robust", {"o": 0.05}, {"robust": 0.05}), - ( - "robust wavelet", - {"o": 0.05, "r": 0.01, "l": 1}, - {"robust": 0.05, "alpha": 0.01, "wavelet": True}, - ), - ("one sweep", {"i": 1}, {"maxiter": 1}), -] +def _centred(x: torch.Tensor, axes: tuple[int, ...], inverse: bool = False) -> torch.Tensor: + """numpy's centred unitary transform, which is what the app's samples are in.""" + shifted = torch.fft.ifftshift(x, dim=axes) + moved = (torch.fft.ifftn if inverse else torch.fft.fftn)(shifted, dim=axes, norm="ortho") + return torch.fft.fftshift(moved, dim=axes) -@pytest.mark.parametrize( - "arguments", - [(t, a) for _, t, a in _POCS_CONFIGURATIONS], - ids=[n for n, _, _ in _POCS_CONFIGURATIONS], -) -def test_the_pocsense_app_is_the_tool_to_the_last_bit(arguments): - flags, ours = arguments - kspace, maps = _cartesian() - tool = bt.pocsense(kspace, maps, **{"i": 10, **flags}) - made = apps.pocsense(kspace, maps, **{"maxiter": 10, **ours}) - assert torch.equal(made, tool), f"maximum difference {float((made - tool).abs().max()):.3e}" + +def _range(samples: torch.Tensor, maps: torch.Tensor, axes: tuple[int, ...]) -> torch.Tensor: + """``E E^H`` written out with torch: combine the coils and spread them again.""" + image = (maps.conj() * _centred(samples, axes, inverse=True)).sum(0, keepdim=True) + return _centred(maps * image, axes) -def test_an_odd_grid_is_the_tool_too(): - """The modulation is a phase rather than a sign on an odd axis, so which - two of the three factors of the sparsity projection are multiplied first - is visible in the last bits.""" - kspace = bt.phantom(25, coils=COILS, kspace=True) +def _in_range(shape: int | list[int]): + """Coil samples ``E x`` of a phantom, fully sampled, and the maps they were made with.""" + kspace = bt.phantom(shape, coils=COILS, kspace=True) maps = bt.ecalib(kspace, maps=1) - mask = torch.zeros(25, dtype=torch.complex64) - mask[::ACCEL] = 1 - mask[25 // 2 - 2 : 25 // 2 + 2] = 1 - kspace = kspace * mask.reshape(25, 1) + axes = tuple(axis for axis in range(-kspace.ndim + 1, 0) if kspace.shape[axis] > 1) + return _range(kspace, maps, axes), maps, axes - tool = bt.pocsense(kspace, maps, i=10, r=0.01, l=1) - ours = apps.pocsense(kspace, maps, maxiter=10, alpha=0.01, wavelet=True) - assert torch.equal(ours, tool) +@pytest.mark.parametrize("shape", [SIZE, [16, 16, 16]], ids=["plane", "volume"]) +def test_samples_the_coils_could_have_made_are_a_fixed_point(shape): + """Every sample measured and every sample in the range of the maps: both + projections leave it where it is.""" + samples, maps, _ = _in_range(shape) + made = apps.pocsense(samples, maps, maxiter=5) + scale = float(samples.abs().max()) + assert float((made - samples).abs().max()) / scale < 1e-4 -def test_a_volume_is_the_tool_too(): - """Three transformed axes rather than two, which is what the app works - out from the sensitivities rather than assuming.""" - size = 16 - kspace = bt.phantom([size, size, size], coils=COILS, kspace=True) - maps = bt.ecalib(kspace, maps=1) - mask = torch.zeros(size, size, 1, dtype=torch.complex64) - mask[torch.randperm(size, generator=torch.Generator().manual_seed(3))[: size // 2]] = 1 - mask[size // 2 - 1 : size // 2 + 1] = 1 - kspace = kspace * mask - tool = bt.pocsense(kspace, maps, i=5, r=0.01, l=1) - ours = apps.pocsense(kspace, maps, maxiter=5, alpha=0.01, wavelet=True) - assert torch.equal(ours, tool) +@pytest.mark.parametrize( + "arguments", + [{}, {"alpha": 0.1}, {"robust": 0.05}, {"maxiter": 1}], + ids=["plain", "l2", "robust", "one sweep"], +) +def test_the_answer_lies_in_the_range_of_the_coils(arguments): + """Without a sparsity projection the last one in a sweep is onto the + coils, or a scaling of it, so the answer is its own projection.""" + kspace, maps = _cartesian() + made = apps.pocsense(kspace, maps, **{"maxiter": 10, **arguments}) + scale = float(made.abs().max()) + assert float((_range(made, maps, (-2, -1)) - made).abs().max()) / scale < 1e-4 def test_the_sweep_is_the_projections_in_order(): diff --git a/tests/test_cli.py b/tests/test_cli.py index 3890d64b..1fcb0ab1 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -41,8 +41,6 @@ def _dataset(tmp_path, monkeypatch): ["pics", "-r", "0.1", "-i", "20"], ["pics", "-l2", "-r", "0.05", "-i", "20"], ["pics", "-R", "W:3:0:0.01", "-i", "20"], - ["pocsense", "-i", "10"], - ["pocsense", "-i", "10", "-r", "0.01", "-l", "1"], ] @@ -96,9 +94,9 @@ def test_an_unknown_command_is_refused(capsys): def test_help_comes_from_the_catalogue(capsys): """BART answers its own help by calling ``exit``, which in this process would end the interpreter.""" - assert main(["pocsense", "--help"]) == 0 + assert main(["pics", "--help"]) == 0 printed = capsys.readouterr().out - assert "Perform POCSENSE reconstruction" in printed + assert "Parallel-imaging compressed-sensing reconstruction" in printed assert "max. number of iterations" in printed diff --git a/tests/test_cuda.py b/tests/test_cuda.py index bf29f4d6..b4c2ce3b 100644 --- a/tests/test_cuda.py +++ b/tests/test_cuda.py @@ -155,12 +155,11 @@ def test_a_tool_that_is_not_kept_on_the_card_still_answers_on_it(): """Crossing to the host is where a tool runs, not what the caller sees.""" from bartorch._dispatch import _ON_DEVICE - assert "pocsense" not in _ON_DEVICE + assert "fftmod" not in _ON_DEVICE n, coils = 32, 2 - img = bt.phantom([n, n], coils=coils) - ksp = bartorch.fft(img, axes=(-2, -1)).cuda() - maps = (torch.ones(1, coils, 1, n, n, dtype=torch.complex64) / coils**0.5).cuda() + img = bt.phantom([n, n], coils=coils).cuda() - out = bt.pocsense(ksp, maps) + out = bartorch.fftmod(img, (-2, -1)) assert out.device.type == "cuda" - torch.testing.assert_close(out.cpu(), bt.pocsense(ksp.cpu(), maps.cpu()), rtol=1e-4, atol=1e-5) + expected = bartorch.fftmod(img.cpu(), (-2, -1)) + torch.testing.assert_close(out.cpu(), expected, rtol=1e-4, atol=1e-5) diff --git a/tests/test_finufft.py b/tests/test_finufft.py index 7a51d997..59c20126 100644 --- a/tests/test_finufft.py +++ b/tests/test_finufft.py @@ -2067,8 +2067,13 @@ def test_a_library_without_its_modules_transforms_at_the_baseline(tmp_path): from bartorch._lib import library_path - alone = tmp_path / library_path().name + # A repaired wheel finds its vendored libraries through $ORIGIN/../bartorch.libs. + home = library_path().parent + (tmp_path / home.name).mkdir() + alone = tmp_path / home.name / library_path().name shutil.copy(library_path(), alone) + for vendored in home.parent.glob("*.libs"): + (tmp_path / vendored.name).symlink_to(vendored, target_is_directory=True) code = ( "import bartorch.tools as bt; from bartorch import _finufft, linop\n" "t = bt.traj(x=32, y=16, r=True)\n" diff --git a/tests/test_openmp.py b/tests/test_openmp.py index b2c14ad0..6ab92bb3 100644 --- a/tests/test_openmp.py +++ b/tests/test_openmp.py @@ -1,11 +1,9 @@ """One OpenMP runtime in the process, whichever of torch and bartorch comes first. -LLVM's runtime ends the process when a second copy of itself initialises, so -on macOS and Windows the library is linked against the runtime torch carries -(cmake/openmp.cmake). Each case runs in a fresh interpreter, exercises -torch's thread pool and FINUFFT's together, and reads the images the process -has loaded. On Linux the toolchain's runtime is used and GNU's tolerates a -second copy, so there the transform is exercised and the count only reported. +The library carries no runtime of its own: it is bound to the one torch +carries (cmake/openmp.cmake), and a wheel vendors none. Each case runs in a +fresh interpreter, exercises torch's thread pool and FINUFFT's together, and +reads the images the process has loaded. """ import subprocess @@ -99,11 +97,10 @@ def _runtimes(order: str) -> list[str]: def test_threaded_torch_and_threaded_finufft_share_one_openmp_runtime(order): runtimes = _runtimes(order) print(order, runtimes) - if sys.platform in ("darwin", "win32"): - assert len(runtimes) == 1, f"{len(runtimes)} OpenMP runtimes loaded: {runtimes}" - assert "torch" in runtimes[0].replace("\\", "/").lower(), ( - f"the one runtime loaded is not torch's: {runtimes[0]}" - ) + assert len(runtimes) == 1, f"{len(runtimes)} OpenMP runtimes loaded: {runtimes}" + assert "torch" in runtimes[0].replace("\\", "/").lower(), ( + f"the one runtime loaded is not torch's: {runtimes[0]}" + ) def test_loading_modifies_no_installed_file(): diff --git a/tests/test_optional_outputs.py b/tests/test_optional_outputs.py index 2eb43560..95cd5dfe 100644 --- a/tests/test_optional_outputs.py +++ b/tests/test_optional_outputs.py @@ -1,12 +1,13 @@ -"""Commands whose output BART marks optional, and the two that write it anyway. +"""Commands whose output BART marks optional, and the two that need it anyway. -A derived wrapper asks for as many output arrays as BART requires, and nine -commands require none. For seven of them that is right: they print their +A derived wrapper asks for as many output arrays as BART requires, and ten +commands require none. For eight of them that is right: they print their answer and write a file only when given a name, through ``anon_cfl`` or behind -an ``if``. Two create it unguarded -- ``mobafit.c:398`` and ``morphop.c:77`` --- and ``create_cfl`` hands the name to ``io_unlink_if_opened``, which calls -``strcmp`` on it. With no name that is a segmentation fault, which takes the -process rather than raising. +an ``if``. ``morphop`` creates it unguarded (``morphop.c:77``), and +``create_cfl`` hands the name to ``io_unlink_if_opened``, which calls +``strcmp`` on it: with no name that is a segmentation fault, which takes the +process rather than raising. ``raga`` refuses an omitted name unless it is +asked to print a search instead. """ import pytest @@ -68,7 +69,7 @@ def test_the_fit_recovers_a_decay_it_was_given(): def test_a_command_that_prints_its_answer_still_prints_it(): # `estdelay` prints the delays when it is given no output name and writes - # a quadratic fit when it is; passing one so that `mobafit` stops crashing + # a quadratic fit when it is; passing one so that `morphop` stops crashing # must not quietly turn the one into the other. traj = bt.traj(x=16, y=8, r=True) kspace = torch.randn(1, 8, 16, 1, dtype=torch.complex64) @@ -78,7 +79,7 @@ def test_a_command_that_prints_its_answer_still_prints_it(): @pytest.mark.parametrize( - "name", ["bench", "estdelay", "estshift", "ismrmrd", "measure", "roistat", "seq"] + "name", ["bench", "estdelay", "estshift", "ismrmrd", "measure", "roistat", "seq", "stl"] ) def test_the_guarded_commands_are_left_alone(name): # Each of these writes through anon_cfl or behind an if, so BART is happy diff --git a/tests/test_options.py b/tests/test_options.py index 43af61c8..b0da0a2e 100644 --- a/tests/test_options.py +++ b/tests/test_options.py @@ -34,7 +34,6 @@ def test_every_option_of_every_command_reaches_its_own_flag(): ("pics", "ist_last", "--ist_last"), ("pics", "fista_pqr", "--fista_pqr"), ("moba", "scale_data", "--scale_data"), - ("seq", "slice_thickness", "--slice_thickness"), # And the ones it does spell with a hyphen, which must still work. ("pics", "no_toeplitz", "--no-toeplitz"), ("pics", "gpu_gridding", "--gpu-gridding"), diff --git a/tests/test_tools.py b/tests/test_tools.py index f0667747..344025b0 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -318,8 +318,8 @@ def test_the_guard_is_only_for_the_ir_mgre_sequence(): def test_a_command_refuses_a_term_its_parser_does_not_know(): kspace, maps = _pics_data() - with pytest.raises(TypeError, match="sqpics does not take"): - bt.sqpics(kspace, maps, R=priors.TotalGeneralizedVariation((-1, -2), 0.01)) + with pytest.raises(TypeError, match="moba does not take"): + bt.moba(kspace, maps, r=priors.TotalGeneralizedVariation((-1, -2), 0.01)) def test_a_derived_wrapper_is_shaped_like_the_command_line():