Adding options to switch kernel at runtime

DiamonDinoia · DiamonDinoia · commit dbca829a1601 · 2025-12-10T17:37:57.000+01:00
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -83,6 +83,8 @@ if(FINUFFT_USE_CPU)
     add_subdirectory(src)
 endif()
 
+add_subdirectory(src/common)
+
 if(FINUFFT_USE_CUDA)
     include(cuda_setup)
     add_subdirectory(src/cuda)
diff --git a/cmake/utils.cmake b/cmake/utils.cmake
@@ -103,7 +103,7 @@ function(enable_asan target)
 endfunction()
 
 function(finufft_link_test target)
-    target_link_libraries(${target} PRIVATE finufft::finufft)
+    target_link_libraries(${target} PRIVATE finufft::finufft finufft::common)
     if(FINUFFT_USE_DUCC0)
         target_compile_definitions(${target} PRIVATE FINUFFT_USE_DUCC0)
     endif()
diff --git a/include/finufft.fh b/include/finufft.fh
@@ -20,6 +20,7 @@ c     DEPRECATED: spread_kerpad kept for ABI compatibility, ignored by library
       real*8 upsampfac
       integer spread_thread, maxbatchsize, spread_nthr_atomic
       integer spread_max_sp_size
+      integer spread_kernel
       integer fftw_lock_fun, fftw_unlock_fun, fftw_lock_data
 
       end type
diff --git a/include/finufft/spreadinterp.h b/include/finufft/spreadinterp.h
@@ -21,9 +21,9 @@ namespace finufft {
 namespace spreadinterp {
 
 template<typename T>
-FINUFFT_EXPORT_TEST int setup_spreader(finufft_spread_opts &opts, T eps, double upsampfac,
-                                       int kerevalmeth, int debug, int showwarn,
-                                       int spreadinterponly, int dim);
+FINUFFT_EXPORT_TEST int setup_spreader(
+    finufft_spread_opts &opts, T eps, double upsampfac, int kerevalmeth, int debug,
+    int showwarn, int spreadinterponly, int dim, int kernel_type = 0);
 
 int spreadcheck(UBIGINT N1, UBIGINT N2, UBIGINT N3, const finufft_spread_opts &opts);
 template<typename T>
diff --git a/include/finufft_common/kernel.h b/include/finufft_common/kernel.h
@@ -6,6 +6,8 @@
 #include <vector>
 
 #include <finufft_common/constants.h>
+#include <finufft_common/utils.h>
+#include <finufft_spread_opts.h>
 
 namespace finufft::kernel {
 
@@ -56,15 +58,34 @@ template<class T, class F> std::vector<T> fit_monomials(F &&f, int n, T a, T b)
   return c;
 }
 
-template<typename T> T evaluate_kernel(T x, T beta, T c) {
-  /* ES ("exp sqrt") kernel evaluation at single real argument:
-      phi(x) = exp(beta.(sqrt(1 - (2x/n_s)^2) - 1)),    for |x| < nspread/2
-     related to an asymptotic approximation to the Kaiser--Bessel, itself an
-     approximation to prolate spheroidal wavefunction (PSWF) of order 0.
-     This is the "reference implementation", used by eg finufft/onedim_* 2/17/17.
-     Rescaled so max is 1, Barnett 7/21/24
+template<typename T> T evaluate_kernel(T x, T beta, T c, int kernel_type = 0) {
+  /* Kernel evaluation at single real argument.
+     kernel_type == 0 : ES ("exp sqrt") kernel (default)
+       phi_ES(x) = exp(beta*(sqrt(1 - c*x^2) - 1))
+     kernel_type == 1 : Kaiser--Bessel (KB) kernel
+       phi_KB(x) = I_0(beta*sqrt(1 - c*x^2)) / I_0(beta)
+     Note: `std::cyl_bessel_i` from <cmath> is used for I_0.
+     Rescaled so max is 1.
   */
+  if (kernel_type == 1) {
+    // Kaiser--Bessel (normalized by I0(beta)). Use std::cyl_bessel_i from <cmath>.
+    const T inner        = std::sqrt(T(1) - c * x * x);
+    const T arg          = beta * inner;
+    const double i0_arg  = ::finufft::common::cyl_bessel_i(0, static_cast<double>(arg));
+    const double i0_beta = ::finufft::common::cyl_bessel_i(0, static_cast<double>(beta));
+    return static_cast<T>(i0_arg / i0_beta);
+  }
+
+  // default to ES
   return std::exp(beta * (std::sqrt(T(1) - c * x * x) - T(1)));
 }
 
+FINUFFT_EXPORT int compute_kernel_ns(double upsampfac, double tol, int kernel_type,
+                                     const finufft_spread_opts &opts);
+
+FINUFFT_EXPORT void initialize_kernel_params(finufft_spread_opts &opts, double upsampfac,
+                                             double tol, int kernel_type);
+
+FINUFFT_EXPORT double sigma_max_tol(double upsampfac, int kernel_type, int max_ns);
+
 } // namespace finufft::kernel
diff --git a/include/finufft_common/utils.h b/include/finufft_common/utils.h
@@ -13,6 +13,9 @@ namespace common {
 FINUFFT_EXPORT_TEST void gaussquad(int n, double *xgl, double *wgl);
 std::tuple<double, double> leg_eval(int n, double x);
 
+// Series implementation of the modified Bessel function of the first kind I_nu(x)
+double cyl_bessel_i(double nu, double x) noexcept;
+
 // helper to generate the integer sequence in range [Start, End]
 template<int Offset, typename Seq> struct offset_seq;
 
diff --git a/include/finufft_mod.f90 b/include/finufft_mod.f90
@@ -21,6 +21,7 @@ module finufft_mod
    real(kind=C_DOUBLE) :: upsampfac
    integer(kind=C_INT) :: spread_thread, maxbatchsize
    integer(kind=C_INT) :: spread_nthr_atomic, spread_max_sp_size
+   integer(kind=C_INT) :: spread_kernel
    integer(kind=C_SIZE_T) :: fftw_lock_fun, fftw_unlock_fun, fftw_lock_data
    !  really, last should be type(C_PTR) :: etc, but fails to print nicely
 
diff --git a/include/finufft_opts.h b/include/finufft_opts.h
@@ -32,6 +32,8 @@ typedef struct finufft_opts { // defaults see finufft_core.cpp:finufft_default_o
   int spread_nthr_atomic; // if >=0, threads above which spreader OMP critical goes
                           // atomic
   int spread_max_sp_size; // if >0, overrides spreader (dir=1) max subproblem size
+  int spread_function;    // (dev only) 0:DEFAULT, (do not change), there is no guarantee
+                          // what non-zero values do and behaviour can change anytime
                           // sphinx tag (don't remove): @opts_end
 
   // User can provide their own FFTW planner lock functions for thread safety
diff --git a/include/finufft_spread_opts.h b/include/finufft_spread_opts.h
@@ -26,6 +26,8 @@ typedef struct finufft_spread_opts {
   double ES_beta;
   double ES_halfwidth;
   double ES_c;
+  // Kernel selector: 0 = ES (default), 1 = Kaiser--Bessel (KB)
+  int kernel_type; // default 0
 } finufft_spread_opts;
 
 #endif // FINUFFT_SPREAD_OPTS_H
diff --git a/matlab/finufft.mw b/matlab/finufft.mw
@@ -103,6 +103,9 @@ $     }
 $     else if (strcmp(fname[ifield],"spread_max_sp_size") == 0) {
 $       oc->spread_max_sp_size = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
 $     }
+$     else if (strcmp(fname[ifield],"spread_kernel") == 0) {
+$       oc->spread_kernel = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
+$     }
 $     else if (strcmp(fname[ifield],"spreadinterponly") == 0) {
 $       oc->spreadinterponly = (int)round(*mxGetPr(mxGetFieldByNumber(om,idx,ifield)));
 $     }
diff --git a/python/finufft/finufft/_finufft.py b/python/finufft/finufft/_finufft.py
@@ -86,6 +86,7 @@ class FinufftOpts(ctypes.Structure):
                       ('maxbatchsize', c_int),
                       ('spread_nthr_atomic', c_int),
                       ('spread_max_sp_size', c_int),
+                      ('spread_kernel', c_int),
                       ('fftw_lock_fun', c_void_p),
                       ('fftw_unlock_fun', c_void_p),
                       ('fftw_lock_data', c_void_p)]
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -23,7 +23,6 @@ set(FINUFFT_SOURCES
     finufft_core.cpp
     c_interface.cpp
     finufft_utils.cpp
-    utils.cpp
 )
 
 if(FINUFFT_BUILD_FORTRAN)
@@ -49,6 +48,7 @@ if(FINUFFT_USE_DUCC0)
 endif()
 
 target_link_libraries(finufft PRIVATE $<BUILD_INTERFACE:finufft_fftlibs xsimd>)
+target_link_libraries(finufft PRIVATE finufft_common)
 if(FINUFFT_USE_OPENMP)
     target_link_libraries(finufft PRIVATE OpenMP::OpenMP_CXX)
     if(NOT FINUFFT_STATIC_LINKING)
@@ -100,4 +100,5 @@ endif()
 
 set(_targets ${INSTALL_TARGETS})
 list(APPEND _targets finufft)
+list(APPEND _targets finufft_common)
 set(INSTALL_TARGETS "${_targets}" PARENT_SCOPE)
diff --git a/src/common/CMakeLists.txt b/src/common/CMakeLists.txt
@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 3.24)
+
+# Build a small static helper library for CPU-only kernel parameter helpers
+# Source is in ../kernel.cpp as requested
+set(FINUFFT_COMMON_SOURCES kernel.cpp utils.cpp)
+
+add_library(finufft_common STATIC ${FINUFFT_COMMON_SOURCES})
+add_library(finufft::common ALIAS finufft_common)
+
+# The public include directory is the top-level include/
+target_include_directories(
+    finufft_common
+    PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include> $<INSTALL_INTERFACE:include>
+)
+
+# Visibility and compile options consistent with finufft
+if(FINUFFT_SHARED_LINKING)
+    target_compile_definitions(finufft_common PRIVATE FINUFFT_DLL)
+endif()
+
+target_compile_features(finufft_common PRIVATE cxx_std_17)
+
+# If you need -lm etc.
+find_library(MATH_LIBRARY m)
+if(MATH_LIBRARY)
+    target_link_libraries(finufft_common PRIVATE ${MATH_LIBRARY})
+endif()
+
+set_target_properties(finufft_common PROPERTIES CXX_VISIBILITY_PRESET hidden VISIBILITY_INLINES_HIDDEN YES)
diff --git a/src/common/kernel.cpp b/src/common/kernel.cpp
@@ -0,0 +1,79 @@
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <finufft_common/common.h>
+#include <finufft_common/kernel.h>
+#include <finufft_spread_opts.h>
+
+namespace finufft::kernel {
+
+int compute_kernel_ns(double upsampfac, double tol, int kernel_type,
+                      const finufft_spread_opts &opts) {
+  // Note: opts is unused here but kept for API consistency if needed later
+  int ns;
+  if (upsampfac == 2.0) {
+    ns = (int)std::ceil(-std::log10(tol / 10.0));
+  } else {
+    ns = (int)std::ceil(
+        -std::log(tol) / (finufft::common::PI * std::sqrt(1.0 - 1.0 / upsampfac)));
+  }
+  ns = std::max(2, ns);
+  return ns;
+}
+
+void initialize_kernel_params(finufft_spread_opts &opts, double upsampfac, double tol,
+                              int kernel_type) {
+  int ns;
+  // Respect any pre-set opts.nspread (e.g., clipped by caller). If it's <=0 compute it.
+  if (opts.nspread > 0) {
+    ns = opts.nspread;
+  } else {
+    ns = compute_kernel_ns(upsampfac, tol, kernel_type, opts);
+  }
+
+  opts.kernel_type  = kernel_type;
+  opts.nspread      = ns; // ensure opts is populated with the (possibly clipped) ns
+  opts.ES_halfwidth = (double)ns / 2.0;
+  opts.ES_c         = 4.0 / (double)(ns * ns);
+
+  if (kernel_type == 1) {
+    // Kaiser-Bessel (KB)
+    // Formula from Beatty et al. 2005
+    double tmp   = (double)ns * (double)ns / (upsampfac * upsampfac);
+    double term2 = (upsampfac - 0.5) * (upsampfac - 0.5);
+    opts.ES_beta = common::PI * std::sqrt(tmp * term2 - 0.8);
+  } else {
+    // Exponential of Semicircle (ES)
+    double betaoverns = 2.30;
+    if (ns == 2)
+      betaoverns = 2.20;
+    else if (ns == 3)
+      betaoverns = 2.26;
+    else if (ns == 4)
+      betaoverns = 2.38;
+
+    if (upsampfac != 2.0) {
+      double gamma = 0.97;
+      betaoverns   = gamma * common::PI * (1.0 - 1.0 / (2.0 * upsampfac));
+    }
+    opts.ES_beta = betaoverns * (double)ns;
+  }
+
+  if (opts.debug) {
+    const char *kname = (kernel_type == 1) ? "KB" : "ES";
+    printf("setup_spreader: using spread kernel type %d (%s)\n", kernel_type, kname);
+    printf("setup_spreader eps=%.3g sigma=%.3g (%s): chose ns=%d beta=%.3g\n", tol,
+           upsampfac, kname, ns, opts.ES_beta);
+  }
+}
+
+double sigma_max_tol(double upsampfac, int kernel_type, int max_ns) {
+  if (upsampfac == 2.0) {
+    return 10.0 * std::pow(10.0, -(double)max_ns);
+  } else {
+    return std::exp(
+        -(double)max_ns * finufft::common::PI * std::sqrt(1.0 - 1.0 / upsampfac));
+  }
+}
+
+} // namespace finufft::kernel
diff --git a/src/common/utils.cpp b/src/common/utils.cpp
diff --git a/src/cuda/CMakeLists.txt b/src/cuda/CMakeLists.txt
diff --git a/src/finufft_core.cpp b/src/finufft_core.cpp
diff --git a/src/spreadinterp.cpp b/src/spreadinterp.cpp
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
diff --git a/test/accuracy_test.cpp b/test/accuracy_test.cpp
