Defer upsampfact to setpts

This allows leveraging density to choose a better value.

Author: DiamonDinoia

CHANGELOG

@@ -37,6 +37,8 @@ Master (working towards v2.5.0), 11/25/25
   the public API and share a unified STL-based implementation for kersum,
   mass and relative-error checks. This common kernel header will also enable
   reuse by the GPU. (Barbone, Barnett, Lu, PR #748)
+* `makeplan` now defers upsampfact choice in `setpts` if upsampfact is set to auto.
+  This allows the decision to account for the `density`.
 
 V 2.4.1 7/8/25
 

docs/opts.rst

@@ -199,15 +199,6 @@ for only two settings, as follows. Otherwise, setting it to zero chooses a good
 * ``upsampfac=1.25`` : low-upsampling option, with lower RAM, smaller FFTs, but wider spreading kernel. The latter can be much faster than the standard when the number of nonuniform points is similar or smaller to the number of modes, and/or if low accuracy is required. It is especially much (2 to 3 times) faster for type 3 transforms. However, the kernel widths :math:`w` are about 50% larger in each dimension, which can lead to slower spreading (it can also be faster due to the smaller size of the fine grid). Because the kernel width is limited to 16, currently, thus only 9-digit accuracy can currently be reached when using ``upsampfac=1.25``.
 
 
-**hint_nj**: Estimated number of nonuniform points available at plan time.
-If nonzero, ``makeplan`` uses this estimate to choose ``upsampfac``.
-Each ``setpts`` call recomputes the factor using the actual ``nj`` and
-re-initializes the plan if it changes.  A value of ``0`` defers the choice
-until ``setpts``.  If ``upsampfac`` is set explicitly, ``hint_nj`` is ignored.
-
-
-
-
 **spread_thread**: in the case of multiple transforms per call (``ntr>1``, or the "many" interfaces), controls how multithreading is used to spread/interpolate each batch of data.
 
 * ``spread_thread=0`` : makes an automatic choice between the below. Recommended.

include/finufft.fh

@@ -16,7 +16,6 @@ c     alg performance opts...
       integer nthreads, fftw, spread_sort, spread_kerevalmeth
       integer spread_kerpad
       real*8 upsampfac
-      integer*4 hint_nj
       integer spread_thread, maxbatchsize, spread_nthr_atomic
       integer spread_max_sp_size
       integer fftw_lock_fun, fftw_unlock_fun, fftw_lock_data

include/finufft/fft.h

@@ -80,6 +80,8 @@ template<> struct Finufft_FFT_plan<float> {
     for (auto i : dims) nf *= i;
 #ifdef _OPENMP
     lock();
+    fftwf_plan_with_nthreads(nthreads);
+#endif
     // If a previous plan exists, destroy it before creating a new one.
     if (plan_) {
       fftwf_destroy_plan(plan_);
@@ -89,12 +91,6 @@ template<> struct Finufft_FFT_plan<float> {
       fftwf_destroy_plan(plan_adj_);
       plan_adj_ = nullptr;
     }
-#else
-    lock();
-#endif
-#ifdef _OPENMP
-    fftwf_plan_with_nthreads(nthreads);
-#endif
     plan_     = fftwf_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
                                     reinterpret_cast<fftwf_complex *>(ptr), nullptr, 1,
                                     int(nf), reinterpret_cast<fftwf_complex *>(ptr), nullptr,
@@ -169,7 +165,9 @@ template<> struct Finufft_FFT_plan<double> {
     for (auto i : dims) nf *= i;
 #ifdef _OPENMP
     lock();
-    // If a previous plan exists, destroy it before creating a new one.
+    fftw_plan_with_nthreads(nthreads);
+#endif
+    // If a previous plan exists, destroy it before creating a new on
     if (plan_) {
       fftw_destroy_plan(plan_);
       plan_ = nullptr;
@@ -178,12 +176,6 @@ template<> struct Finufft_FFT_plan<double> {
       fftw_destroy_plan(plan_adj_);
       plan_adj_ = nullptr;
     }
-#else
-    lock();
-#endif
-#ifdef _OPENMP
-    fftw_plan_with_nthreads(nthreads);
-#endif
     plan_     = fftw_plan_many_dft(int(dims.size()), dims.data(), int(batchSize),
                                    reinterpret_cast<fftw_complex *>(ptr), nullptr, 1, int(nf),
                                    reinterpret_cast<fftw_complex *>(ptr), nullptr, 1, int(nf),

include/finufft_mod.f90

@@ -17,7 +17,6 @@ module finufft_mod
    integer(kind=C_INT) :: nthreads,fftw,spread_sort,spread_kerevalmeth
    integer(kind=C_INT) :: spread_kerpad
    real(kind=C_DOUBLE) :: upsampfac
-   integer(kind=C_INT) :: hint_nj
    integer(kind=C_INT) :: spread_thread, maxbatchsize
    integer(kind=C_INT) :: spread_nthr_atomic, spread_max_sp_size
    integer(kind=C_SIZE_T) :: fftw_lock_fun, fftw_unlock_fun, fftw_lock_data

include/finufft_opts.h

@@ -5,8 +5,6 @@
 #ifndef FINUFFT_OPTS_H
 #define FINUFFT_OPTS_H
 
-#include <stddef.h>
-
 typedef struct finufft_opts { // defaults see finufft_core.cpp:finufft_default_opts_t()
   // sphinx tag (don't remove): @opts_start
   // FINUFFT options:
@@ -28,7 +26,6 @@ typedef struct finufft_opts { // defaults see finufft_core.cpp:finufft_default_o
   int spread_kerevalmeth; // spreader: 0 exp(sqrt()), 1 Horner piecewise poly (faster)
   int spread_kerpad;      // (exp(sqrt()) only): 0 don't pad kernel to 4n, 1 do
   double upsampfac;       // upsampling ratio sigma: 2.0 std, 1.25 small FFT, 0.0 auto
-  int hint_nj;            // estimated nj at plan time; 0 means unknown
   int spread_thread;      // (vectorized ntr>1 only): 0 auto, 1 seq multithreaded,
                           //                          2 parallel single-thread spread
   int maxbatchsize;       // (vectorized ntr>1 only): max transform batch, 0 auto

python/finufft/finufft/_finufft.py

@@ -13,7 +13,6 @@
 from ctypes import c_int
 from ctypes import c_longlong
 from ctypes import c_void_p
-from ctypes import c_size_t
 
 import numpy as np
 import os
@@ -83,7 +82,6 @@ class FinufftOpts(ctypes.Structure):
                       ('spread_kerevalmeth', c_int),
                       ('spread_kerpad', c_int),
                       ('upsampfac', c_double),
-                      ('hint_nj', c_int),
                       ('spread_thread', c_int),
                       ('maxbatchsize', c_int),
                       ('spread_nthr_atomic', c_int),

src/finufft_core.cpp

@@ -588,7 +588,6 @@ void finufft_default_opts_t(finufft_opts *o)
   o->spread_kerevalmeth = 1;
   o->spread_kerpad      = 1;
   o->upsampfac          = 0.0;
-  o->hint_nj            = 0;
   o->spread_thread      = 0;
   o->maxbatchsize       = 0;
   o->spread_nthr_atomic = -1;
@@ -788,22 +787,9 @@ FINUFFT_PLAN_T<TF>::FINUFFT_PLAN_T(int type_, int dim_, const BIGINT *n_modes, i
 
   if (tol < std::numeric_limits<TF>::epsilon()) ier = FINUFFT_WARN_EPS_TOO_SMALL;
 
-  // If user hasn't forced upsampfac and hint_nj > 0 (types 1,2 only), choose an
-  // initial upsampfac now using the hint. The actual upsampfac may still change
-  // at setpts if the real nj differs from the hint. If hint_nj == 0, defer
-  // upsampfac selection to setpts. For type 3, hint is ignored (always defers).
-  if (!upsamp_locked && type != 3 && opts.hint_nj > 0) {
-    double density = double(opts.hint_nj) / double(N());
-    opts.upsampfac = bestUpsamplingFactor<TF>(opts.nthreads, density, dim, type, tol);
-    if (opts.debug > 1)
-      printf("[%s] chosen upsampfac=%.2f from hint (density=%.3g)\n", __func__,
-             opts.upsampfac, density);
-    // Note: don't set upsamp_locked - the actual setpts may still adjust
-  }
-
-  // If upsampfac is known (forced, or from hint for types 1,2), init spreader/FFT now.
-  // If type 3 with upsampfac=0, or types 1,2 with hint_nj==0, defer to setpts.
-  if (upsamp_locked || opts.upsampfac != 0) {
+  // If user hasn't forced upsampfac, defer selection to setpts.
+  if (upsamp_locked) {
+    // If upsampfac is known (forced), init spreader/FFT now.
     int ier2 = init_spreader_and_fft();
     if (ier2 > 1) throw int(ier2);
     ier = std::max(ier, ier2);
@@ -812,45 +798,6 @@ FINUFFT_PLAN_T<TF>::FINUFFT_PLAN_T(int type_, int dim_, const BIGINT *n_modes, i
   if (ier > 1) throw int(ier); // report setup_spreader status (could be warning)
 }
 
-template<typename TF> void FINUFFT_PLAN_T<TF>::precompute_horner_coeffs() {
-  const auto nspread   = spopts.nspread;
-  const auto upsampfac = opts.upsampfac;
-  const auto max_degree =
-      upsampfac == 2.00 ? horner_nc_200(nspread) : horner_nc_125(nspread);
-  // get the xsimd padding
-  static constexpr auto simd_size = xsimd::batch<TF>::size;
-  const auto padded_ns            = (nspread + simd_size - 1) & -simd_size;
-  // pad degree if desired; keep equal to max_degree for now
-  const auto padded_degree = max_degree;
-  horner_coeffs.fill(0);
-
-  for (int j = 0; j < nspread; ++j) {
-    const auto kernel = [&](const TF x) {
-      // scale manually from -1, 1 to -w/2+j to -w/2+j+1
-      return evaluate_kernel(TF(0.5) * TF(x - nspread + 2 * j + 1),
-                             TF(this->spopts.ES_beta), TF(this->spopts.ES_c));
-    };
-    // disable polyfit scaling
-    static constexpr TF a = -1.0;
-    static constexpr TF b = 1.0;
-    const auto coeffs     = fit_monomials(kernel, max_degree, a, b);
-    // store transposed: for each degree k
-    for (size_t k = 0; k < coeffs.size(); ++k)
-      horner_coeffs[k * padded_ns + j] = coeffs[k];
-  }
-
-  if (opts.debug > 2) {
-    // Print transposed layout: all degree 0 coeffs for intervals, then degree 1, ...
-    for (size_t k = 0; k < padded_degree; ++k) {
-      printf("[%s] degree=%lu: ", __func__, k);
-      for (size_t j = 0; j < padded_ns; ++j) // use padded_ns to show padding as well
-        printf("%g ", horner_coeffs[k * padded_ns + j]);
-      printf("\n");
-    }
-  }
-}
-
-// MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
 template<typename TF>
 int finufft_makeplan_t(int type, int dim, const BIGINT *n_modes, int iflag, int ntrans,
                        TF tol, FINUFFT_PLAN_T<TF> **pp, const finufft_opts *opts)

test/CMakeLists.txt

@@ -11,7 +11,6 @@ set(TESTS
     finufft3dkernel_test
     spreadinterp1d_test
     adjointness
-    hint_nj_test
 )
 
 foreach(TEST ${TESTS})
@@ -99,8 +98,6 @@ function(add_tests_with_prec PREC REQ_TOL CHECK_TOL SUFFIX)
     )
 
     add_test(NAME run_adjointness_${PREC} COMMAND adjointness${SUFFIX} WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
-
-    add_test(NAME run_hint_nj_${PREC} COMMAND hint_nj_test${SUFFIX} WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
 endfunction()
 
 # use above function to actually add the tests, with certain requested and check