Skip to content

Improve min sigma bounds#845

Open
DiamonDinoia wants to merge 4 commits into
flatironinstitute:masterfrom
DiamonDinoia:ups_bounds_estimator1
Open

Improve min sigma bounds#845
DiamonDinoia wants to merge 4 commits into
flatironinstitute:masterfrom
DiamonDinoia:ups_bounds_estimator1

Conversation

@DiamonDinoia
Copy link
Copy Markdown
Collaborator

@DiamonDinoia DiamonDinoia commented Apr 3, 2026

Since we now can now use any sigma we should also account for rounding error in to prevent giving wrong results. I started with @mreineck suggestion in #841 (comment) then refined it using emprical data. Asked Claude to cleanup all my scratch code into seomthing sensible. The final plots are below. It is not perfect but I think is good enough for the purpose. When N is small next235 influences sigma indirectly so any fitting breaks. With larger N the fit improves and follows the data more closely.
sigma_fit_double
sigma_fit_float

Edit: updated the plots with the original @mreineck formula.

@DiamonDinoia DiamonDinoia requested a review from lu1and10 April 3, 2026 20:04
@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 3, 2026

That's super interesting! I was aware that my model was crude, but I had expected it to be too pessimistic rather than too optimistic - I need to do some more digging.

@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 3, 2026

Ah, OK, if you don't use the N-dependent floor, then the results are somewhat understandable. But without that term, the model doesn't make much sense in for the scenarios you are plotting...

@DiamonDinoia
Copy link
Copy Markdown
Collaborator Author

DiamonDinoia commented Apr 3, 2026

Ah, OK, if you don't use the N-dependent floor, then the results are somewhat understandable. But without that term, the model doesn't make much sense in for the scenarios you are plotting...

Hang on: const double gridlen = *std::max_element(mstu.begin(), mstu.begin() + dim); and eps_round = 0.48 * eps_mach * gridlen accounts for that no?

@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 3, 2026

I was looking at _mreineck_sigma_min in find_sigma_bound.py and didn't see this kind of line, just the comment "mreineck's 3-term model: kernel + rdyn (no N-dependent floor)." And that seems to match with the graphs.

@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 3, 2026

If you look at the float32, type1 graphs, the estimated sigmas from the "mreineck" model all reach 2.0 between 1e-7 and 1e-6, independent of N. With the N-dependent term, they should arrive there at very different epsilons.

@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 3, 2026

Also, I think you need to use the maximum element of nf123, not mstu.

DiamonDinoia added a commit to DiamonDinoia/finufft that referenced this pull request Apr 4, 2026
Address PR flatironinstitute#845 review: use fine grid size (nfdim) instead of mode
count (mstu) for the rounding floor in the sigma estimator. Move
check from makeplan to setpts where nfdim is known. Restore
mreineck's original 3-term formula (with phase error term) in
the comparison plots.
@DiamonDinoia DiamonDinoia force-pushed the ups_bounds_estimator1 branch from 555a756 to b65ad95 Compare April 4, 2026 03:39
@DiamonDinoia
Copy link
Copy Markdown
Collaborator Author

If you look at the float32, type1 graphs, the estimated sigmas from the "mreineck" model all reach 2.0 between 1e-7 and 1e-6, independent of N. With the N-dependent term, they should arrive there at very different epsilons.

This is Claude showing his limitations. Among the sea of code I produced to fit the polynomial it forgot which was actually the original formula. I plotted the 3 contributions separately to understand what was going on by zeroing one of them a the time and left the code with the last one as 0.

@DiamonDinoia
Copy link
Copy Markdown
Collaborator Author

Also, I think you need to use the maximum element of nf123, not mstu.

Also, the check should have been done in setpts as I do not know N in makeplan. Refining the existing warning threw me off.

@mreineck
Copy link
Copy Markdown
Collaborator

mreineck commented Apr 4, 2026

Thank you, that looks much closer to what I would expect!

Comment thread src/common/kernel.cpp
return std::min(1.0 / (1.0 - u * u), MAXSIGMA);
}

double estimated_tol(double sigma, int ns, int dim, double eps_mach, double gridlen) {
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this function used anywhere, or is it designed for later usage?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is not used anymore. I did use it at some point. Arguably not necessary now. But, I think it can be useful to further refine this estimator in the future.

Copy link
Copy Markdown
Member

@lu1and10 lu1and10 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good. Do we want to support the ES kernel in this PR?

DiamonDinoia added a commit to DiamonDinoia/finufft that referenced this pull request Apr 22, 2026
Address PR flatironinstitute#845 review: use fine grid size (nfdim) instead of mode
count (mstu) for the rounding floor in the sigma estimator. Move
check from makeplan to setpts where nfdim is known. Restore
mreineck's original 3-term formula (with phase error term) in
the comparison plots.
@DiamonDinoia DiamonDinoia force-pushed the ups_bounds_estimator1 branch from b65ad95 to 8eb66ce Compare April 22, 2026 16:27
DiamonDinoia and others added 3 commits April 22, 2026 12:56
Replace the old polynomial-based SigmaEstimator with a calibrated
two-term analytical model for predicting the minimum upsampling factor
(sigma) that achieves a requested tolerance.

The model: tol = tolfac * exp(-(ns-0.1)*pi*u^1.19) + 117*eps_mach
where u = sqrt(1-1/sigma). Fits empirical sigma_min to within 0.03
sigma for double precision, dim 1, types 1-3.

- Add find_sigma_bound.cpp: C++ training tool that binary-searches
  for empirical sigma_min at each tolerance
- Add find_sigma_bound.py: Python validation/calibration script with
  PR comparison plots (analytical, mreineck 3-term, calibrated)
- Update makeplan.hpp: new estimated_tol using calibrated formula,
  lowest_sigma_impl returns MAXSIGMA (not MAXSIGMA+1) when not
  achievable, warning always suggests a concrete sigma value

Co-authored-by: Artur <magisterbrownie@gmail.com>
Hybrid kernel-aliasing / rounding-floor model predicts the minimum
upsampling factor needed for a given tolerance. Check moved to setpts
where nfdim (fine grid size) is known, using it as gridlen instead of
mstu. Throws FINUFFT_ERR_EPS_TOO_SMALL (or warns if allow_eps_too_small).

New functions: estimated_tol, lowest_sigma in finufft::common.
Includes Python validation tool devel/find_sigma_bound.py.
Address PR flatironinstitute#845 review: use fine grid size (nfdim) instead of mode
count (mstu) for the rounding floor in the sigma estimator. Move
check from makeplan to setpts where nfdim is known. Restore
mreineck's original 3-term formula (with phase error term) in
the comparison plots.
@DiamonDinoia DiamonDinoia force-pushed the ups_bounds_estimator1 branch 2 times, most recently from c21477e to 3d67408 Compare April 22, 2026 18:32
check_sigma previously force-threw FINUFFT_ERR_EPS_TOO_SMALL in the
"rounding floor dominates" branch even when the user had explicitly
opted in via opts.allow_eps_too_small=1. That negated the opt-in's
contract.

Library change: gate the throw solely on !allow_eps_too_small. The
warning still distinguishes the unachievable case (no upsampfac can
help) from the curable case (print the suggested upsampfac floor).

Test changes kept minimal:
- adjointness.cpp: set allow_eps_too_small=1 once at the top.
- basicpassfail.cpp: loosen tol to 1e-3 for single prec only, since
  the single-prec rounding floor at N=1e3 already exceeds 1e-5.
- test/CMakeLists.txt: keep float tol at 1e-3 (pending auto-sigma work).

All 26 ctests pass (float + double).
@DiamonDinoia DiamonDinoia force-pushed the ups_bounds_estimator1 branch from 3d67408 to b4045e5 Compare April 22, 2026 20:21
@ahbarnett ahbarnett self-assigned this May 5, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants