|
| 1 | +/* test/tolsweep: pass-fail accuracy test for either float/double CPU FINUFFT |
| 2 | + that sweeps across full range of tolerances, dims, types, for set of upsampfacs |
| 3 | + (just two std USFs for now). |
| 4 | + Uses relative L2 error norms, with direct reference evaluation. |
| 5 | + Exit code: zero if success, nonzero upon failure. |
| 6 | +
|
| 7 | + Based on Barbone's accuracy_test and Barnett's matlab/test/tolsweeptest.m |
| 8 | + The logic is taken from the latter (no significance to decades). Barnett 1/5/26 |
| 9 | + Multiple dims, USFs, extra debug output for failures. 1/6/26. |
| 10 | + nthr=1 (for valgrind!) and tweak for 6s (float), 12s (double) runtimes. 1/7/26. |
| 11 | +*/ |
| 12 | + |
| 13 | +#include <cmath> |
| 14 | +#include <cstdlib> |
| 15 | +#include <finufft/test_defs.h> |
| 16 | +#include <iomanip> |
| 17 | +#include <iostream> |
| 18 | +#include <vector> |
| 19 | +// test utilities: direct DFT and norm helpers |
| 20 | +#include "utils/dirft1d.hpp" |
| 21 | +#include "utils/dirft2d.hpp" |
| 22 | +#include "utils/dirft3d.hpp" |
| 23 | +#include "utils/norms.hpp" |
| 24 | + |
| 25 | +int main(int argc, char *argv[]) { |
| 26 | + |
| 27 | + // Define test problems, tolerance ranges, slack factors... |
| 28 | + BIGINT M = 500; // pick problem size: # sources (balance runtime vs rand-averaging) |
| 29 | + // N vectors to test: first triplet is for dim=1, then for dim=2, etc... |
| 30 | + BIGINT Nm_alldims[3][3] = {{30, 1, 1}, {20, 40, 1}, {10, 11, 12}}; // Ntot~1e3 ok |
| 31 | + int ntr = 1; // >1 in tolsweeptest.m but only for speed/convenience |
| 32 | + int isign = +1; |
| 33 | + |
| 34 | + // *** make these vary by dim?: had to grow type-3 slack=15 here for macos CI... |
| 35 | + double tolslack[3] = {5.0, 5.0, 15.0}; // tunable slack parameters for each type |
| 36 | + double tolsperdecade = 8; // controls overall effort (tol resolution) |
| 37 | + double tolstep = pow(10.0, -1.0 / tolsperdecade); // multiplicative tol step, <1 |
| 38 | + constexpr FLT EPSILON = std::numeric_limits<FLT>::epsilon(); |
| 39 | + double mintol = 0.5 * EPSILON; // somewhat arbitrary where start (catch warns) |
| 40 | + int ntols = std::ceil(log(mintol) / log(tolstep)); |
| 41 | + |
| 42 | + // Defaults |
| 43 | + int kerformula = 0; |
| 44 | + int showwarn = 0; |
| 45 | + int verbose = 2; |
| 46 | + int debug = 0; |
| 47 | + // test set of upsampfacs each with matching error floor for each dim... |
| 48 | + const int nu = 2; // how many USFs (just the standard ones for now) |
| 49 | + double upsampfac[nu] = {1.25, 2.0}; |
| 50 | +#ifdef SINGLE |
| 51 | + double floor[nu][3] = {{1e-4, 1e-3, 1e-2}, {1e-5, 1e-5, 1e-5}}; // inner is dim |
| 52 | +#else |
| 53 | + double floor[nu][3] = {{3e-9, 3e-9, 1e-8}, {3e-14, 3e-14, 3e-14}}; |
| 54 | +#endif |
| 55 | + |
| 56 | + // If user asked for help, print usage and exit |
| 57 | + for (int ai = 1; ai < argc; ++ai) { |
| 58 | + if (std::string(argv[ai]) == "-h" || std::string(argv[ai]) == "--help") { |
| 59 | + std::cout << "Usage: " << argv[0] << " [kerformula [showwarn [verbose [debug]]]]\n"; |
| 60 | + std::cout |
| 61 | + << " kerformula : spread kernel formula (0:default, >0: for experts)\n"; |
| 62 | + std::cout |
| 63 | + << " showwarn : whether to print warnings (0=silent default, 1=show)\n"; |
| 64 | + std::cout |
| 65 | + << " verbose : 0 silent, 1 summary, 2 +fail deets (default), 3...\n"; |
| 66 | + std::cout << " debug : passed to opts.debug for every call\n"; |
| 67 | + std::cout << "Example: " << argv[0] << " 0 1 2 0\n"; |
| 68 | + return 0; |
| 69 | + } |
| 70 | + } |
| 71 | + if (argc > 1) sscanf(argv[1], "%d", &kerformula); |
| 72 | + if (argc > 2) sscanf(argv[2], "%d", &showwarn); |
| 73 | + if (argc > 3) sscanf(argv[3], "%d", &verbose); |
| 74 | + if (argc > 4) sscanf(argv[4], "%d", &debug); |
| 75 | + |
| 76 | + finufft_opts opts{}; |
| 77 | + FINUFFT_DEFAULT_OPTS(&opts); |
| 78 | + opts.spread_kerformula = kerformula; |
| 79 | + opts.debug = debug; |
| 80 | + opts.showwarn = showwarn; |
| 81 | + opts.nthreads = 1; // single-threaded FINUFFT faster since small (esp valgrind!) |
| 82 | + |
| 83 | + std::vector<FLT> x(M), y(M), z(M), X, Y, Z; // xyz real vs XYZ freq-space |
| 84 | + std::vector<CPX> c(M), ce(M), F, Fe; // (we don't know N yet, since varies by dim) |
| 85 | + srand(42); // fix seed |
| 86 | + int nfailtot = 0; // overall count across all dims, USF, tols, types |
| 87 | + |
| 88 | + for (int dim = 1; dim <= 3; ++dim) { /////////////////////// loop over dims |
| 89 | + if (verbose) printf("%s: %dD =============================\n", argv[0], dim); |
| 90 | + BIGINT *Nm = Nm_alldims[dim - 1]; // ptr to Nm array (3-el) for this dim |
| 91 | + BIGINT N = Nm[0] * Nm[1] * Nm[2]; // tot # modes, or freq-pts for type 3 |
| 92 | + X.resize(N); |
| 93 | + Y.resize(N); |
| 94 | + Z.resize(N); |
| 95 | + F.resize(N); |
| 96 | + Fe.resize(N); |
| 97 | + |
| 98 | + for (int u = 0; u < nu; ++u) { // ========================= loop over upsampfacs |
| 99 | + // (sigma) |
| 100 | + opts.upsampfac = upsampfac[u]; |
| 101 | + if (verbose) printf(" upsampfac = %.3g -----------------\n", opts.upsampfac); |
| 102 | + |
| 103 | + double worstfac[3] = {0}; // track the largest clearance for each type |
| 104 | + double tol = 1.0; // starting (max) tol to test |
| 105 | + int npass[3] = {0}, nfail[3] = {0}; // counts for each type |
| 106 | + for (int t = 0; t < ntols; ++t) { // ............... loop over tols |
| 107 | + |
| 108 | + for (int type = 1; type <= 3; ++type) { // ---------------------- loop over types |
| 109 | + |
| 110 | + // just make new data each test & type, even data not needed for that type |
| 111 | + for (BIGINT j = 0; j < M; ++j) { |
| 112 | + x[j] = PI * randm11(); |
| 113 | + y[j] = PI * randm11(); |
| 114 | + z[j] = PI * randm11(); |
| 115 | + c[j] = crandm11(); |
| 116 | + } |
| 117 | + for (BIGINT k = 0; k < N; ++k) { |
| 118 | + X[k] = Nm[0] * rand01(); // type 3: scale freq s,t,u NU pts by "mode" sizes |
| 119 | + Y[k] = Nm[1] * rand01(); // (should match erralltypedim.m) |
| 120 | + Z[k] = Nm[2] * rand01(); |
| 121 | + F[k] = crandm11(); |
| 122 | + } |
| 123 | + FINUFFT_PLAN plan; // do tested transform... |
| 124 | + int ier = FINUFFT_MAKEPLAN(type, dim, Nm, isign, ntr, (FLT)tol, &plan, &opts); |
| 125 | + int ier_set = FINUFFT_SETPTS(plan, M, x.data(), y.data(), z.data(), N, X.data(), |
| 126 | + Y.data(), Z.data()); |
| 127 | + FINUFFT_EXECUTE(plan, c.data(), F.data()); // type 2 writes to c, others to F |
| 128 | + FINUFFT_DESTROY(plan); |
| 129 | + |
| 130 | + if (dim == 1) // do the relevant (of nine) direct "exact" evals... |
| 131 | + if (type == 1) |
| 132 | + dirft1d1<BIGINT>(M, x, c, isign, Nm[0], Fe); // exact ans written into Fe |
| 133 | + else if (type == 2) |
| 134 | + dirft1d2<BIGINT>(M, x, ce, isign, Nm[0], F); // exact ans written into ce |
| 135 | + else |
| 136 | + dirft1d3<BIGINT>(M, x, c, isign, Nm[0], X, Fe); // exact ans written into Fe |
| 137 | + else if (dim == 2) |
| 138 | + if (type == 1) |
| 139 | + dirft2d1<BIGINT>(M, x, y, c, isign, Nm[0], Nm[1], Fe); |
| 140 | + else if (type == 2) |
| 141 | + dirft2d2<BIGINT>(M, x, y, ce, isign, Nm[0], Nm[1], F); |
| 142 | + else |
| 143 | + dirft2d3<BIGINT>(M, x, y, c, isign, N, X, Y, Fe); |
| 144 | + else // dim=3 |
| 145 | + if (type == 1) |
| 146 | + dirft3d1<BIGINT>(M, x, y, z, c, isign, Nm[0], Nm[1], Nm[2], Fe); |
| 147 | + else if (type == 2) |
| 148 | + dirft3d2<BIGINT>(M, x, y, z, ce, isign, Nm[0], Nm[1], Nm[2], F); |
| 149 | + else |
| 150 | + dirft3d3<BIGINT>(M, x, y, z, c, isign, N, X, Y, Z, Fe); |
| 151 | + |
| 152 | + double relerr; // compute relevant error metric |
| 153 | + if (type == 2) |
| 154 | + relerr = relerrtwonorm<BIGINT>(M, ce, c); // ||ce-c||/||ce|| so ce comes 1st |
| 155 | + else |
| 156 | + relerr = relerrtwonorm<BIGINT>(N, Fe, F); |
| 157 | + |
| 158 | + if (ier > 1 || ier_set > 1) { // an error, not merely warning, we exit |
| 159 | + fprintf(stderr, " tolsweep: %dD%d failed! ier=%d, ier_setpts=%d\n", dim, |
| 160 | + type, ier, ier_set); |
| 161 | + return 1; |
| 162 | + } |
| 163 | + |
| 164 | + if (ier == 0) { |
| 165 | + int ti = type - 1; // index for 3-el arrays |
| 166 | + double req = std::max(floor[u][dim - 1], tolslack[ti] * tol); // threshold |
| 167 | + double clearfac = relerr / req; // factor by which beats req (<=1 ok, >1 fail) |
| 168 | + worstfac[ti] = std::max(worstfac[ti], clearfac); // track the worst case |
| 169 | + bool pass = (relerr <= req); |
| 170 | + if (pass) { |
| 171 | + ++npass[ti]; |
| 172 | + if (verbose > 2) |
| 173 | + printf( |
| 174 | + " %dd%d, tol %8.3g:\trelerr = %.3g, \tclearancefac=%.3g\tpass\n", |
| 175 | + dim, type, tol, relerr, clearfac); |
| 176 | + } else { |
| 177 | + ++nfail[ti]; |
| 178 | + printf( |
| 179 | + " %dd%d, tol %8.3g:\trelerr = %.3g, \tclearancefac=%.3g \tFAIL\n", |
| 180 | + dim, type, tol, relerr, clearfac); |
| 181 | + if (verbose > 1) { |
| 182 | + printf(" Rerunning with " |
| 183 | + "debug=1:_______________________________________\n"); |
| 184 | + opts.debug = 1; |
| 185 | + FINUFFT_MAKEPLAN(type, dim, Nm, isign, ntr, (FLT)tol, &plan, &opts); |
| 186 | + FINUFFT_SETPTS(plan, M, x.data(), y.data(), z.data(), N, X.data(), |
| 187 | + Y.data(), Z.data()); |
| 188 | + FINUFFT_EXECUTE(plan, c.data(), F.data()); // type 2 writes to c |
| 189 | + FINUFFT_DESTROY(plan); |
| 190 | + printf(" (Rerun " |
| 191 | + "done)__________________________________________________\n"); |
| 192 | + opts.debug = debug; // reset to cmdline arg value |
| 193 | + } |
| 194 | + } |
| 195 | + } else // finufft returned warning (ie cannot achieve accuracy): don't test |
| 196 | + if (verbose > 2) |
| 197 | + printf( |
| 198 | + " %dd%d, tol %8.3g:\trelerr = %.3g, \t(warn ier=%d: not tested)\n", |
| 199 | + dim, type, tol, relerr, ier); |
| 200 | + |
| 201 | + } // --------------------------- |
| 202 | + |
| 203 | + tol *= tolstep; // reduce tol in geometric progr |
| 204 | + } // ........................................... |
| 205 | + |
| 206 | + if (verbose) |
| 207 | + for (int ti = 0; ti < 3; ++ti) // cout beats printf buffer-flush (for |
| 208 | + // check_finufft.sh) |
| 209 | + std::cout << " " << dim << "D" << ti + 1 << " summary: " << npass[ti] |
| 210 | + << " pass, " << nfail[ti] << " fail. worstfac=" << worstfac[ti] |
| 211 | + << std::endl; |
| 212 | + nfailtot += nfail[0] + nfail[1] + nfail[2]; |
| 213 | + |
| 214 | + } // ========================== |
| 215 | + |
| 216 | + } ////////////////////////// |
| 217 | + |
| 218 | + return nfailtot > 0; // if any fails, allows all test cases to complete |
| 219 | +} |
0 commit comments