Claude/jitfields fastfields migration v5r416

MIGRATION.md

@@ -26,17 +26,27 @@ Each ported operation appears at three "library" levels:
 | module | kernels | cpu-impl | cuda-impl | cpu-lib | cuda-lib | lib | CPU tested |
 |----------------|:------:|:--------:|:---------:|:-------:|:--------:|:---:|:----------:|
 | distance | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | **yes** |
-| posdef | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
+| posdef | ✓ | ✓ | ✓ | ✓ | ~ | ✓ | **yes** |
+| resize | ✓ | ✓ | ✓ | ✓ | ~ | ✓ | **yes** |
+| restrict | ✓ | ✓ | ✓ | ✓ | ~ | ✓ | **yes** |
+| splinc | ✓ | ✓ | ✓ | ✓ | ~ | ✓ | **yes** |
 | pushpull | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
 | reg_field | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
 | reg_flow | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
-| resize | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
-| restrict | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
-| splinc | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no |
 | tetrahedron | ✓ | ✓ | — | ✗ | ✗ | ✗ | no |
 
-"✓" for a lib column means the dispatch layer exists; only `distance` is currently
-compiled and tested (CPU). `tetrahedron` has no cuda-impl header yet.
+"✓" for a cpu-lib/lib column means the dispatch layer exists **and is CPU-compiled+tested**.
+"~" for cuda-lib means the dispatch source is written and mirrors cpu-lib, but is **not yet
+compilable**: the corresponding `cuda-impl/<module>.h` only defines device (`CUGLOB`) kernels
+and lacks the host launchers (allocate/copy shape+stride to device, launch, forward `stream`)
+that the cuda-lib layer calls — see the "not yet fixed" section. Such modules are intentionally
+**left out of the cuda-lib Makefile `MODULES`** until the launchers exist, so the CUDA build
+keeps working. `tetrahedron` has no cuda-impl header yet.
+
+### Public API naming
+`resize`/`restrict`/`splinc` would be an illegal same-name namespace+function inside `ff::cpu`,
+so the public ops are named **`resample`** (resize), **`restriction`** (restrict), and
+**`spline_coeff`** (splinc). `restriction` accumulates into `out`, so callers pre-zero it.
 
 ## Build & test (CPU)
 
@@ -64,9 +74,37 @@ CPU path against a brute-force / reference implementation, as `test_distance.cpp
  the L1 transform was never called.
 5. **cpu-lib/Makefile, lib/Makefile** — object compile rule missing `-fPIC`, so the
  shared-library link failed (`relocation R_X86_64_PC32 … recompile with -fPIC`).
+6. **kernels/threadpool.inl** — defined non-`inline` free functions and two
+ namespace-scope globals (`internal::num_threads`, `internal::global_pool`) in a
+ header. Fine for a single-module library, but a library with ≥2 module objects
+ (e.g. distance + posdef) failed to link (`multiple definition of …`). Made the
+ helpers `inline` and wrapped the two globals in Meyers-singleton accessors (C++11
+ has no inline variables). This unblocks every multi-module library.
+7. **kernels/posdef/utils.h** — three alias templates referenced undeclared names
+ (`_as_points`, `left`/`right`, and `_return_type` with no args); any posdef include
+ failed to parse. **kernels/posdef/sym.inl** — dynamic `Sym::invert` called the
+ inherited `copy_` unqualified (two-phase lookup) → `this_type::copy_`.
+8. **cpu-impl/posdef.h** — runtime `nbatch` passed as a template arg to `index2offset`;
+ a dead static-C `else` branch that misbound the dynamic specialization under C++11
+ two-arm instantiation (guarded with `(C<0?1:C)`); `delete`→`delete[]` at 4 sites.
+9. **kernels/utils.h** — `prod<size>(x)` called `typed_prod<T,size>(x, size)` with a
+ spurious extra arg (no matching overload); triggered by `restrict::loop`.
+10. **cpu-impl/{resize,restrict,splinc}.h** — wrong include prefix `"lib/…"` →
+ `"kernels/…"`; impl namespace was plain `ff::<module>` but the kernels live in
+ `ff::cpu::` (`FF_DEVICE`) so it must be `FF_NAMESPACE_BEGIN(FF)/(FF_DEVICE)/(<module>)`
+ like distance; `index2offset_nd<ndim>()` runtime-ndim → dynamic overload;
+ `jf::has_atomic_add` → `has_atomic_add`.
 
 ## Bugs found, NOT yet fixed (need a CUDA build to verify)
 
+- **cuda-impl/{posdef,resize,restrict,splinc}.h** — provide only device (`CUGLOB`)
+ kernels and **no host launchers** analogous to distance's `CUHOST dt()`. The
+ cuda-lib layer for these modules is written but cannot compile until host wrappers
+ are added, so these modules are omitted from the cuda-lib Makefile `MODULES` for
+ now. Additionally cuda-impl uses the same wrong `"lib/…"` includes, and
+ `resize.h`'s `kernelnd` passes an undefined `x` (should be `loc`) to
+ `Multiscale<ndim>::resize`. (Tasks tracked separately.)
+
 - **cuda-impl/distance_euclidean.h** — `dt()` allocates a scratch buffer of
  `vector_size * (sizeof(offset_t) + 2*sizeof(scalar_t))` bytes, but `dt_kernel`
  addresses `z`/`d` at `buf + stride_buf * n * …` and indexes each of the

Makefile

@@ -24,7 +24,7 @@ MOSUF = o
 SOSUF = so
 SONAME = soname
 OMPFLAG = -fopenmp
-RPATH = -Wl,-rpath,$$ORIGIN/../lib
+RPATH = -Wl,-rpath,'$$ORIGIN'/../lib
 USE_OPENMP ?= 0
 
 ########################################################################
@@ -108,7 +108,13 @@ clean-cpu:
 
 MODULES = \
 	distance \
-	posdef
+	posdef \
+	resize \
+	restrict \
+	splinc \
+	pushpull \
+	reg_field \
+	reg_flow
 
 OBJECTS = $(addprefix $(BUILDDIR)/,$(addsuffix .$(MOSUF),$(MODULES)))
 CPPFILES = $(addsuffix .cpp,$(MODULES))

distance.cpp

@@ -5,7 +5,7 @@
 #include "cuda/distance.h"
 #endif
 
-using namespace FF;
+FF_NAMESPACE_BEGIN(FF)
 
 #define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
 #define IS_CPU(tensor) (tensor.device.device_type == DLDeviceType::kDLCPU || \
@@ -123,3 +123,5 @@ void dt_mesh(
 
  throw std::invalid_argument("unsupported device");
 }
+
+FF_NAMESPACE_END(FF)

distance.h

@@ -5,6 +5,8 @@
 
 FF_NAMESPACE_BEGIN(FF)
 
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
 FF_NAMESPACE_BEGIN(bound_t)
 using T = int8_t;
 static constexpr T Dynamic = -1; ///< Used to turn-off static implementations in templated classes
@@ -30,6 +32,7 @@ static constexpr T FifthOrder = 5;
 static constexpr T SixthOrder = 6;
 static constexpr T SeventhOrder = 7;
 FF_NAMESPACE_END(spline_t)
+#endif // FF_LIB_BOUND_SPLINE_T
 
 /**
  * @brief Compute the Euclidean distance transform of a tensor.

posdef.cpp

@@ -5,7 +5,7 @@
 #include "cuda/posdef.h"
 #endif
 
-using namespace FF;
+FF_NAMESPACE_BEGIN(FF)
 
 #define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
 #define IS_CPU(tensor) (tensor.device.device_type == DLDeviceType::kDLCPU || \
@@ -136,3 +136,5 @@ void sym_invert_(
 
  throw std::invalid_argument("unsupported device");
 }
+
+FF_NAMESPACE_END(FF)

pushpull.cpp

@@ -0,0 +1,90 @@
+#include <stdexcept>
+#include "pushpull.h"
+#include "cpu/pushpull.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/pushpull.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor) (tensor.device.device_type == DLDeviceType::kDLCPU || \
+ tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void pull(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline,
+ int8_t bound,
+ int8_t extrapolate,
+ int stream )
+{
+#ifdef FF_WITH_CUDA
+ if (IS_CUDA(out))
+ return FF_CUDA::pull(out, inp, grid, spline, bound, extrapolate, stream);
+#endif
+ if (IS_CPU(out))
+ return FF_CPU::pull(out, inp, grid, spline, bound, extrapolate, stream);
+
+ throw std::invalid_argument("unsupported device");
+}
+
+void push(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline,
+ int8_t bound,
+ int8_t extrapolate,
+ int stream )
+{
+#ifdef FF_WITH_CUDA
+ if (IS_CUDA(out))
+ return FF_CUDA::push(out, inp, grid, spline, bound, extrapolate, stream);
+#endif
+ if (IS_CPU(out))
+ return FF_CPU::push(out, inp, grid, spline, bound, extrapolate, stream);
+
+ throw std::invalid_argument("unsupported device");
+}
+
+void count(
+ DLTensor & out,
+ const DLTensor & grid,
+ int8_t spline,
+ int8_t bound,
+ int8_t extrapolate,
+ int stream )
+{
+#ifdef FF_WITH_CUDA
+ if (IS_CUDA(out))
+ return FF_CUDA::count(out, grid, spline, bound, extrapolate, stream);
+#endif
+ if (IS_CPU(out))
+ return FF_CPU::count(out, grid, spline, bound, extrapolate, stream);
+
+ throw std::invalid_argument("unsupported device");
+}
+
+void grad(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline,
+ int8_t bound,
+ int8_t extrapolate,
+ bool abs,
+ int stream )
+{
+#ifdef FF_WITH_CUDA
+ if (IS_CUDA(out))
+ return FF_CUDA::grad(out, inp, grid, spline, bound, extrapolate, abs, stream);
+#endif
+ if (IS_CPU(out))
+ return FF_CPU::grad(out, inp, grid, spline, bound, extrapolate, abs, stream);
+
+ throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)

pushpull.h

@@ -0,0 +1,118 @@
+#ifndef FF_LIB_PUSHPULL
+#define FF_LIB_PUSHPULL
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic = -1;
+static constexpr T Zero = 0;
+static constexpr T Replicate = 1;
+static constexpr T DCT1 = 2;
+static constexpr T DCT2 = 3;
+static constexpr T DST1 = 4;
+static constexpr T DST2 = 5;
+static constexpr T DFT = 6;
+static constexpr T NoCheck = 7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic = -1;
+static constexpr T Nearest = 0;
+static constexpr T Linear = 1;
+static constexpr T Quadratic = 2;
+static constexpr T Cubic = 3;
+static constexpr T FourthOrder = 4;
+static constexpr T FifthOrder = 5;
+static constexpr T SixthOrder = 6;
+static constexpr T SeventhOrder = 7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Sample ("pull") a spline-encoded volume at arbitrary coordinates.
+ *
+ * Channel-last, x-first coordinate convention:
+ * inp : (*batch, *inshape, C)
+ * grid : (*batch, *outshape, D) with D == the spatial rank (1, 2 or 3)
+ * out : (*batch, *outshape, C)
+ *
+ * @param out Output tensor (pulled samples)
+ * @param inp Input volume (spline coefficients)
+ * @param grid Sampling coordinates (in voxels, x-first)
+ * @param spline Spline order applied to every spatial dim
+ * @param bound Boundary condition applied to every spatial dim
+ * @param extrapolate 1: always; 0: not past voxel centres; -1: not past edges
+ * @param stream Cuda stream on which to operate
+ */
+void pull(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline = spline_t::Quadratic,
+ int8_t bound = bound_t::DCT2,
+ int8_t extrapolate = 1,
+ int stream = 0
+);
+
+/**
+ * @brief Splat ("push") values into a volume; numerical adjoint of `pull`.
+ * `out` must be pre-zeroed by the caller (values are accumulated).
+ *
+ * inp : (*batch, *outshape, C)
+ * grid : (*batch, *outshape, D)
+ * out : (*batch, *inshape, C)
+ */
+void push(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline = spline_t::Quadratic,
+ int8_t bound = bound_t::DCT2,
+ int8_t extrapolate = 1,
+ int stream = 0
+);
+
+/**
+ * @brief Splat ones (== push of an all-ones input). `out` (*batch,*inshape,1)
+ * must be pre-zeroed.
+ */
+void count(
+ DLTensor & out,
+ const DLTensor & grid,
+ int8_t spline = spline_t::Quadratic,
+ int8_t bound = bound_t::DCT2,
+ int8_t extrapolate = 1,
+ int stream = 0
+);
+
+/**
+ * @brief Sample the spatial gradients of a spline-encoded volume.
+ *
+ * inp : (*batch, *inshape, C)
+ * grid : (*batch, *outshape, D)
+ * out : (*batch, *outshape, C, D)
+ */
+void grad(
+ DLTensor & out,
+ const DLTensor & inp,
+ const DLTensor & grid,
+ int8_t spline = spline_t::Quadratic,
+ int8_t bound = bound_t::DCT2,
+ int8_t extrapolate = 1,
+ bool abs = false,
+ int stream = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_PUSHPULL