From 9078c647f484061bf7d1737dfcd83e0680ae5f8a Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 07:59:55 +0000 Subject: [PATCH 1/7] Move BART to Codeberg master with the rebased downstream stack The submodule moves from 694ab5fe on bartorch-v1.0.00 to 7fab96d7 on downstream: Codeberg master 3beb5572 plus the 22-commit portability stack. What the library needs from that BART: - C23, which BART now compiles as (bool is a keyword, stdbool.h is gone). - USE_GPU beside USE_CUDA, the macro BART's device code is now under, and simu/gpu_bloch.cu in the CUDA sources. - nufft_create2 and sense_nc_init take a field map and a time map. Both are only served by BART's explicit transform (--nufft-conf dft), which is handed straight to BART; a field map without it is refused by name. sense.c includes modelnc.h so the compiler holds the substitute to the header. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- CMakeLists.txt | 6 +++--- external/bart | 2 +- src/csrc/abi/api.c | 2 ++ src/csrc/ops/ops.c | 2 +- src/csrc/ops/sense.c | 17 ++++++++++++----- src/csrc/substitute/nufft_finufft.c | 21 ++++++++++++++++----- src/csrc/substitute/psf.c | 4 ++-- 7 files changed, 37 insertions(+), 17 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7533f908..488dabec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ 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) @@ -171,7 +171,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 +269,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() 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/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]; From b3efb6a13505c1a5191105b45169d4b4d04f6b5d Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 07:59:55 +0000 Subject: [PATCH 2/7] Regenerate the catalogue for the new BART and place its commands BART no longer has pocsense, sqpics or rtnlinv, so their derived tools go, with the CLI's pocsense adapter and the tests that ran them. apps.pocsense stays, built as before from BART's pocs iteration, and is held to a fully sampled phantom being its fixed point and to its answer lying in the range of the coils. New commands are private: pulseq reads and writes files, sort is an array operation torch has, bet and extractdc are not wrapped yet. unwrap takes a set of axes rather than one; mobafit now requires its output and raga refuses to run without one, so the optional-output exception moves from mobafit to raga. The catalogue reader learns OPTL_FLVEC7 and a help string declared as a pointer. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- AGENTS.md | 26 +- docs/api/apps.md | 2 +- docs/api/tools.md | 3 - docs/explanation/non-cartesian.md | 4 +- scripts/gen_catalogue.py | 4 +- src/bartorch/_call.py | 20 +- src/bartorch/_catalogue.py | 517 ++++++++++++++---------------- src/bartorch/_coverage.py | 7 +- src/bartorch/_dispatch.py | 2 +- src/bartorch/apps/pocsense.py | 8 +- src/bartorch/cli/_apps.py | 22 -- src/bartorch/tools/recon.py | 3 - src/bartorch/util.py | 2 +- tests/test_apps.py | 84 ++--- tests/test_cli.py | 6 +- tests/test_cuda.py | 11 +- tests/test_optional_outputs.py | 19 +- tests/test_options.py | 1 - tests/test_tools.py | 4 +- 19 files changed, 329 insertions(+), 416 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index e7860e6c..27dacf8b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -314,7 +314,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 +612,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 +893,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 +926,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 +1101,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 +1121,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/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/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/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..373ce372 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: 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/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/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_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(): From a6a370e173fac5cd7d6d61b6bebdfaa8db999f35 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 08:05:03 +0000 Subject: [PATCH 3/7] Escape the pipes in BART's option help mobafit's -P help writes |M|, which reStructuredText reads as a substitution, and the docs build treats that warning as an error. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- src/bartorch/_call.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/bartorch/_call.py b/src/bartorch/_call.py index 373ce372..d1e3af2e 100644 --- a/src/bartorch/_call.py +++ b/src/bartorch/_call.py @@ -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."] From 87e42582ef72d832a3c212b21557818c0b38d84d Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 08:13:36 +0000 Subject: [PATCH 4/7] Give the copied library its vendored libraries in the module-less FINUFFT test Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- tests/test_finufft.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) 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" From 1801694d49d556e8bf39744ee7408adf77217bc9 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 08:40:35 +0000 Subject: [PATCH 5/7] Keep BART's BARTLIB_API symbols out of the dynamic table The new misc/dllspec.h gives BARTLIB_API default visibility outside Windows, which put seq_*, num_rand_init and debug_printf in the library's dynamic table beside the bartorch_* ABI. Defining the header's guard with BARTLIB_API empty leaves them hidden, as Windows already has them through BARTLIB_STATIC, and a test holds the table's C names to the ABI. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- CMakeLists.txt | 4 ++++ tests/test_abi.py | 23 +++++++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 488dabec..f2644d63 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -467,6 +467,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 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 == [] From 40866e733a34ab6079fed5d87d071f17d1e44bfe Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 08:57:10 +0000 Subject: [PATCH 6/7] Build the CUDA sources with hidden visibility as well nvcc's host code took no visibility preset, so the CUDA build exported BART's cuda_* wrappers, eigenmapscu and wl3_cuda_* beside the bartorch_* ABI, which the dynamic-table test now catches. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index f2644d63..1bb5909e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,6 +33,7 @@ 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) From 58b06cc0fdc405c94f0deb4321df69997c4077c0 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 09:31:05 +0000 Subject: [PATCH 7/7] Bind the Linux library to torch's libgomp and vendor none The Linux wheels carried a libgomp of their own, so a process that imported torch had two OpenMP runtimes. torch carries libgomp.so.1 in torch/lib from 2.7.1 on, which is the name the library's NEEDED entry asks for: the library is loaded after torch on Linux as it already is on macOS and Windows, the wheels are repaired with --exclude libgomp.so.1, and $ORIGIN/../torch/lib is on the runpath. Linux requires torch 2.7.1, the first release whose copy has that name, and the OpenMP test now holds every platform to one runtime, torch's. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014ok1jumatvU8YSzRCEaj8Q --- .github/workflows/publish.yml | 17 +++++++++++------ AGENTS.md | 15 ++++++++++----- CMakeLists.txt | 7 +++++++ cmake/openmp.cmake | 4 +++- docs/design/finufft-embedding.md | 10 ++++++---- docs/guides/user/installation.md | 13 +++++++------ docs/guides/user/prerequisites.md | 2 +- docs/misc/license.md | 1 - pyproject.toml | 3 +++ src/bartorch/_lib.py | 11 +++++------ tests/test_openmp.py | 19 ++++++++----------- 11 files changed, 61 insertions(+), 41 deletions(-) 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 27dacf8b..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.** diff --git a/CMakeLists.txt b/CMakeLists.txt index 1bb5909e..741787ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -511,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/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/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/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/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/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():