diff --git a/Project.toml b/Project.toml index b3d26f2..1d7f8e2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "JetPack" uuid = "24ef3835-3876-54c3-8a7a-956cf69ca0b2" -version = "2.2.1" +version = "2.3.0" [deps] FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" diff --git a/docs/JopRBF/JopRBF-demo.md b/docs/JopRBF/JopRBF-demo.md new file mode 100644 index 0000000..31110d1 --- /dev/null +++ b/docs/JopRBF/JopRBF-demo.md @@ -0,0 +1,158 @@ +# JopRBF: scattered-node model parameterization for JetPack + +`JopRBF` maps a small set of scattered control-node coefficients to a regular fine +model grid using **normalized compact-support radial basis functions** (Wendland +C2). It is a reduced (coarse) model parameterization for `Jets` inverse problems +(FWI dimension reduction, regularization, preconditioning), with an exact adjoint +and a matrix-free implementation that scales to very large grids (for example +1000^3). + +This directory holds the operator's demo scripts and the figures below. The demos +use a local project (`docs/JopRBF/Project.toml`) that adds PyPlot alongside a +dev'd JetPack, so the plotting dependency stays out of JetPack itself: + +``` +julia --startup-file=no --project=docs/JopRBF -t 8 docs/JopRBF/JopRBF_waterbottom_demo.jl # water-bottom comparison figure +julia --startup-file=no --project=docs/JopRBF -t 8 docs/JopRBF/JopRBF_demo.jl # kernel / smoothness figures +``` + +## The operator + +$$ +d(x) = \frac{\sum_j \varphi\!\left(r_j(x)\right)\, c_j}{\sum_j \varphi\!\left(r_j(x)\right)}, +\qquad +r_j(x) = \sqrt{\sum_k \left(\frac{x_k - \xi_{j,k}}{\delta_{k,j}}\right)^2} +$$ + +$$ +\varphi(r) = (1-r)^4\,(4r+1)\ \text{ for } 0 \le r < 1, \qquad \varphi(r) = 0 \text{ otherwise (Wendland C2).} +$$ + +Here $c$ are the node coefficients (the optimization variables), $\xi_j$ the node +locations (in fine-grid index coordinates), and $\delta$ the compact-support +radius. It may be a scalar, a per-axis vector (for anisotropic grids), or a +per-node $\delta_{k,j}$ matrix, so each node can carry its own support radius. +Properties, each a unit test in `test/jop_RBF.jl`: + +- **Compact support, so it scales.** $\varphi$ vanishes for $r \ge 1$, so each + fine point only sees the $O(1)$ nodes within $\delta$. The forward and adjoint + are applied matrix-free with a node bucket index (no dense $N_\text{fine}\times M$ + matrix is formed). The forward threads over fine points and gathers nearby + nodes; the adjoint threads over nodes and walks each node's fine-grid box. Both + enumerate the identical set of (fine point, node) pairs, so the adjoint is exact. +- **C2 smooth, no hot zones.** $\varphi$, $\varphi'$ and $\varphi''$ all vanish at + $r = 1$, so the field is twice continuously differentiable (no bilinear-style + creases). +- **Reproduces constants.** The Shepard normalization is a partition of unity. +- **No overshoot.** $\varphi \ge 0$, so the field is a convex combination of + nearby coefficients and stays within their min and max (no ringing / halos). +- **Build once, reuse.** The node bucket index and the normalization field + $1/\sum_j\varphi$ depend only on the node positions and $\delta$; they are + computed once at construction and reused for every apply. There is no linear + solve. + +### API + +```julia +# nodes are a (ndims(range) x M) matrix in fine-grid index coordinates +A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, nz, nx), nodes; delta = 20.0) +d = A * rand(domain(A)) # coarse coefficients -> fine model +g = A' * rand(range(A)) # adjoint -> gather to nodes + +# anisotropic support (per axis) +A = JopRBF(JetSpace(Float32, M), JetSpace(Float32, nz, nx, ny), nodes3; delta = [30.0, 50.0, 50.0]) + +# per-node support radii (ndims x M): small support shallow, large support deep +A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, nz, nx), nodes; delta = deltas) +``` + +## Water-bottom "freeze" workflow + +The scattered nodes make it natural to parameterize only the sediment below a +(possibly dipping) water bottom while freezing the water column. + +1. Get the water-bottom depth `wb(x)` (for example from `water_bottom_index`). +2. Place the RBF centers below the water bottom with a depth-tapered scattered + cloud ([`tapered_rbf_nodes`](rbf_tapered_nodes.jl), pure base Julia, no external + mesher). The rows **conform to** (run parallel to) the water bottom: each row + sits at a fixed offset from the per-column WB and the offsets march downward by + the local spacing $h(z) = \lambda(z)/\text{ppw}(z)$, $\lambda(z) = v(z)/f$, so + the node spacing follows a points-per-wavelength rule (denser where the + wavelength is short), with `ppw` tapering linearly with depth (dense at the + water bottom, sparse at depth). Making the rows parallel to the WB avoids + horizontal rows cross-cutting a dipping boundary. One capped "ghost" row sits + just **above** the water bottom (masked out of the model, but its support + reaches below) so the near-WB fit has two-sided support and does not overshoot + at the boundary. +3. Build the reduced parameterization operator as a composition: + + ```julia + A = JopRBF(JetSpace(T, M), JetSpace(T, nz, nx), nodes; delta) + S = JopDiagonal(Float64.(below)) # 1 below the water bottom, 0 in the water + P = S ∘ A # reduced parameterization (a JopLn) + ``` + +4. The full model is the frozen water plus the RBF sediment: + $v(c) = v_\text{water} \cdot (\text{water mask}) + P\,c$. + + In FWI you invert for $c$. The reduced gradient is $A'$ applied to the masked + fine-grid gradient ($P' = A' S'$), so the water column is never updated: the + freeze is automatic because there are no nodes in the water and `S` zeros the + water rows. To start from a given sediment model, fit the coefficients once, + `c = P \ vec(v_start .* below)`, then invert. Bounds and regularization now act + on the few hundred coefficients `c` instead of every fine sample, and the + sediment is smooth (C2) and free of overshoot by construction. + +**Figure 1**, produced by +[`JopRBF_waterbottom_demo.jl`](JopRBF_waterbottom_demo.jl), shows this on a 101×201 +model with a sloping water bottom, comparing three inversion frequencies (lower +frequency $\Rightarrow$ longer wavelength $\Rightarrow$ coarser node cloud). Row 1: +the true model, and the partition-of-unity coverage below the water bottom (water +shown white). Rows 2-4: the frozen-water + RBF parameterization at 6, 3 and 1.5 Hz +(nodes overlaid) next to the signed relative error $\Delta v/v$ (red/blue, shared +color scale). Node counts follow the points-per-wavelength taper: 6 Hz ≈ 608 nodes +(25× reduction, 0.0008 RMS), 3 Hz ≈ 186 (83×, 0.0013), 1.5 Hz ≈ 65 (238×, 0.0032) +below-WB relative error. Lower frequency = fewer nodes = a smoother, coarser +representation. + +![Figure 1: water bottom freeze workflow](images/jopRBF_waterbottom_2d.png) + +## Kernel and smoothness + +**Figure 2** shows a single normalized basis function (a compact bump), the +Wendland C2 kernel next to a C0 linear "tent" kernel, and a slice confirming the +basis is exactly zero beyond $\pm\delta$. + +![Figure 2: basis and kernel](images/jopRBF_basis_2d.png) + +**Figure 3** explains why the C2 kernel matters. Both panels interpolate the *same* +coefficients on the *same* coarse node cloud (white dots), differing only in the +kernel. The linear "tent" kernel gives a **C0** field (continuous value but +discontinuous slope) with visible creases/kinks between nodes; the Wendland kernel +gives a smooth **C2** field (continuous value, slope, and curvature), which is what +`JopRBF` uses. The transect along the red line (right) makes it explicit: the C0 +profile kinks at node boundaries, the C2 profile is smooth. Those C0 kinks are the +"hot zones" a piecewise-linear (bilinear or barycentric) parameterization produces +between nodes. + +![Figure 3: C0 vs C2 smoothness](images/jopRBF_smoothness_2d.png) + +## Why this design + +- **Global RBF** (thin-plate, multiquadric) has dense interpolation matrices and + global support, so it does not scale to 1000^3. +- **Triangulation plus C1 cubic** (Clough-Tocher) is local and smooth but needs a + 3D Delaunay mesh and fragile C1 macro-elements. +- **Barycentric-linear on a mesh** is cheap but only C0, which reintroduces the + creases / hot zones. + +Normalized compact-support RBF keeps locality (scales), is C2 (no hot zones), is +meshless (no 3D triangulation), reproduces constants, and does not overshoot. + +## Tests + +`test/jop_RBF.jl`: Wendland C2 kernel (value and first/second derivatives vanish +at the support edge, non-negative); dot-product (adjoint) tests in 1D/2D/3D; +partition of unity; compact support (influence vanishes beyond $\delta$); no +overshoot (convex combination); anisotropic per-axis $\delta$; per-node +(multiresolution) $\delta$; Float32. diff --git a/docs/JopRBF/JopRBF_demo.jl b/docs/JopRBF/JopRBF_demo.jl new file mode 100644 index 0000000..eb4d08c --- /dev/null +++ b/docs/JopRBF/JopRBF_demo.jl @@ -0,0 +1,143 @@ +# Kernel / smoothness figures for JopRBF (normalized compact-support Wendland-C2 RBF). +# +# Generates two PNGs under docs/JopRBF/images/: +# 1. jopRBF_basis_2d.png a single normalized basis, its C2 vs C0 kernel, compact support +# 2. jopRBF_smoothness_2d.png why C2: linear (C0) interpolation kinks between nodes; Wendland (C2) is smooth +# +# Run (global environment has JetPack dev'd + PyPlot): +# julia --startup-file=no docs/JopRBF/JopRBF_demo.jl +# +ENV["MPLBACKEND"] = "Agg" +using JetPack, Jets, LinearAlgebra, Random, Printf, PyPlot + +Random.seed!(20260709) +const OUT = joinpath(@__DIR__, "images") +isdir(OUT) || mkpath(OUT) + +# scattered node cloud (denser shallow) used for the basis figure +function make_nodes(nz, nx) + pts = Vector{Float64}[] + for x in range(1, nx, length = 12) + push!(pts, [1.0, x]); push!(pts, [Float64(nz), x]) + end + for z in range(1, nz, length = 9) + push!(pts, [z, 1.0]); push!(pts, [z, Float64(nx)]) + end + for _ = 1:120 + push!(pts, [1 + (nz - 1) * rand()^0.65, 1 + (nx - 1) * rand()]) + end + reduce(hcat, pts) +end + +# normalized (Shepard) RBF evaluation on the full grid, for a chosen kernel +function eval_normalized(nodes, c, nz, nx, delta, kernel) + d = zeros(nz, nx) + M = size(nodes, 2) + for iz = 1:nz, ix = 1:nx + num = 0.0; den = 0.0 + for j = 1:M + r = sqrt((iz - nodes[1, j])^2 + (ix - nodes[2, j])^2) / delta + w = kernel(r) + num += w * c[j]; den += w + end + d[iz, ix] = den > 0 ? num / den : 0.0 + end + d +end + +wend_c2(r) = r < 1 ? (1 - r)^4 * (4r + 1) : 0.0 # Wendland C2 kernel (what JopRBF uses) +linear0(r) = r < 1 ? (1 - r) : 0.0 # linear "tent" kernel -> C0 (piecewise linear) + +# --------------------------------------------------------------------------- +# Figure 1: a single normalized basis + C2 vs C0 kernel + compact support +# --------------------------------------------------------------------------- +nz, nx = 160, 220 +delta = 26.0 +nodes = make_nodes(nz, nx) +M = size(nodes, 2) +A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, nz, nx), nodes; delta = delta) + +jc = argmin(vec(sum((nodes .- [80.0, 110.0]) .^ 2, dims = 1))) +e = zeros(domain(A)); e[jc] = 1.0 +basis = A * e +zc, xc = nodes[1, jc], nodes[2, jc] + +figure(figsize = (13, 4.2)) +subplot(1, 3, 1) +im = imshow(basis, aspect = "auto", cmap = "viridis") +scatter([xc], [zc], s = 45, c = "red", marker = "x") +xlim(1, nx); ylim(nz, 1) +title("normalized basis A e_j (compact bump)"); xlabel("x"); ylabel("z"); colorbar(im, fraction = 0.046) + +subplot(1, 3, 2) +rr = range(0, 1.3, length = 400) +plot(rr, wend_c2.(rr), lw = 2.4, label = "Wendland C2 (1-r)^4(4r+1)") +plot(rr, linear0.(rr), lw = 1.8, ls = "--", label = "linear C0 (1-r)") +axvline(1.0, color = "k", lw = 0.7, alpha = 0.5) +xlabel("r = |x - ξ| / δ"); ylabel("φ(r)"); title("kernel shape: C2 vs C0"); legend(); grid(alpha = 0.3) + +subplot(1, 3, 3) +row = clamp(round(Int, zc), 1, nz) +plot(1:nx, basis[row, :], lw = 2.0) +axvline(xc - delta, color = "r", ls = ":", lw = 1); axvline(xc + delta, color = "r", ls = ":", lw = 1) +xlabel("x"); ylabel("basis value"); title("slice at z=$(row): zero beyond ±δ"); grid(alpha = 0.3) +tight_layout() +savefig(joinpath(OUT, "jopRBF_basis_2d.png"), dpi = 120); close() + +# --------------------------------------------------------------------------- +# Figure 2: why C2. Same coarse nodes and coefficients, two interpolation +# kernels. Linear weights give a piecewise-linear (C0) field that kinks between +# nodes; the Wendland kernel (what JopRBF uses) is smooth (C2). +# --------------------------------------------------------------------------- +gz, gx = 130, 180 +h = 26.0 +snodes = let pts = Vector{Float64}[] + zs = collect(range(1.0, gz, step = h)); (gz - zs[end] > 1) && push!(zs, Float64(gz)) + xs = collect(range(1.0, gx, step = h)); (gx - xs[end] > 1) && push!(xs, Float64(gx)) + for z in zs, x in xs + jz = (1 < z < gz) ? 0.18h * (rand() - 0.5) : 0.0 + jx = (1 < x < gx) ? 0.18h * (rand() - 0.5) : 0.0 + push!(pts, [z + jz, x + jx]) + end + reduce(hcat, pts) +end +sM = size(snodes, 2) +sdelta = 1.15h +crand = rand(sM) +dC0 = eval_normalized(snodes, crand, gz, gx, sdelta, linear0) +dC2 = eval_normalized(snodes, crand, gz, gx, sdelta, wend_c2) +cmin = min(minimum(dC0), minimum(dC2)); cmax = max(maximum(dC0), maximum(dC2)) +row = 66 + +figure(figsize = (15, 4.4)) +subplot(1, 3, 1) +im = imshow(dC0, aspect = "auto", cmap = "viridis", vmin = cmin, vmax = cmax) +scatter(snodes[2, :], snodes[1, :], s = 12, c = "white", edgecolors = "k", linewidths = 0.3) +axhline(row, color = "r", lw = 1.5) +xlim(1, gx); ylim(gz, 1) +title("linear kernel → C0 field\n(kinks/creases between nodes)"); xlabel("x"); ylabel("z"); colorbar(im, fraction = 0.046) + +subplot(1, 3, 2) +im = imshow(dC2, aspect = "auto", cmap = "viridis", vmin = cmin, vmax = cmax) +scatter(snodes[2, :], snodes[1, :], s = 12, c = "white", edgecolors = "k", linewidths = 0.3) +axhline(row, color = "r", lw = 1.5) +xlim(1, gx); ylim(gz, 1) +title("Wendland C2 kernel → C2 field\n(smooth; what JopRBF uses)"); xlabel("x"); ylabel("z"); colorbar(im, fraction = 0.046) + +subplot(1, 3, 3) +plot(1:gx, dC0[row, :], lw = 1.8, ls = "--", color = "tab:orange", label = "C0 linear (kinks)") +plot(1:gx, dC2[row, :], lw = 2.6, color = "tab:blue", label = "C2 Wendland (smooth)") +xlabel("x"); ylabel("field value along red line"); title("transect at z=$(row)") +legend(); grid(alpha = 0.3) +tight_layout() +savefig(joinpath(OUT, "jopRBF_smoothness_2d.png"), dpi = 120); close() + +# scaling note +let x = rand(domain(A)), y = rand(range(A)) + A * x; A' * y + t1 = @elapsed for _ = 1:5; A * x; end + t2 = @elapsed for _ = 1:5; A' * y; end + @printf("timing on %dx%d, M=%d, %d threads: forward %.1f ms, adjoint %.1f ms\n", + nz, nx, M, Threads.nthreads(), 1000t1 / 5, 1000t2 / 5) +end +println("wrote jopRBF_basis_2d.png and jopRBF_smoothness_2d.png in ", OUT) diff --git a/docs/JopRBF/JopRBF_waterbottom_demo.jl b/docs/JopRBF/JopRBF_waterbottom_demo.jl new file mode 100644 index 0000000..182c8c9 --- /dev/null +++ b/docs/JopRBF/JopRBF_waterbottom_demo.jl @@ -0,0 +1,131 @@ +#= +Figure generation for JopRBF in a water-bottom "freeze" workflow. + +Builds a 101 x 201 (nz x nx) velocity model with a sloping water bottom, parameterizes only the sediment +below the water bottom with a depth-tapered scattered RBF node cloud, and freezes the water column. Produces +docs/JopRBF/images/jopRBF_waterbottom_2d.png. + +Node placement is done by the hand-rolled `tapered_rbf_nodes` (rbf_tapered_nodes.jl, pure base Julia, NO +external mesher): brick-offset rows below the water bottom with spacing(z) = v_trend(z) / (freq * ppw(z)), +ppw linear in depth (dense shallow, sparse deep). The figure compares the parameterization at several +inversion frequencies (lower frequency = longer wavelength = coarser node spacing). Only PyPlot lives in +this directory's local project (docs/JopRBF/Project.toml), not in JetPack itself. + +The reduced parameterization operator is P = S_below ∘ A_RBF , where A_RBF maps node coefficients to the +fine grid and S_below (a JopDiagonal mask) restricts the update to below the water bottom. The full model is + v(c) = v_water · (water mask) + P c +so the water column is frozen (never a function of c). + +Run (uses the local project alongside this script): + julia --startup-file=no --project=docs/JopRBF -t 8 docs/JopRBF/JopRBF_waterbottom_demo.jl +=# +ENV["MPLBACKEND"] = "Agg" +using JetPack, Jets, LinearAlgebra, Random, Printf, PyPlot +include(joinpath(@__DIR__, "rbf_tapered_nodes.jl")) + +Random.seed!(20260709) +const OUT = joinpath(@__DIR__, "images") +isdir(OUT) || mkpath(OUT) + +nz, nx = 101, 201 +vwater = 1500.0 +dz = dx = 10.0 # grid spacing [m] +ppw_top = 6.0 # points per wavelength at the water bottom (dense, shallow) +ppw_bot = 4.0 # points per wavelength at the model bottom (sparse, deep) +support = 1.6 # node support radius = support * local node spacing +node_ms = 22 # node marker size in the figure [points^2] +freqs = (6.0, 3.0, 1.5) # inversion frequencies to compare (lower freq -> coarser node spacing) [Hz] +flabel(f) = @sprintf("%g", f) # frequency label without trailing zeros (6.0 -> "6", 1.5 -> "1.5") + +# sloping water bottom: 12 samples deep on the left, 38 on the right +wbdepth(ix) = 12.0 + (38.0 - 12.0) * (ix - 1) / (nx - 1) + +# "true" model: constant water above the water bottom, v(z) sediment plus two smooth anomalies below it +function true_model() + v = fill(vwater, nz, nx) + for ix = 1:nx, iz = 1:nz + wb = wbcol[ix] # integer WB, CONSISTENT with the below-WB mask (no water leaking into sediment) + if iz >= wb + s = 1650 + 3.2 * (iz - wb) + s += 180 * exp(-(((iz - 55) / 16)^2 + ((ix - 70) / 28)^2)) + s -= 140 * exp(-(((iz - 82) / 18)^2 + ((ix - 150) / 34)^2)) + v[iz, ix] = s + end + end + v +end + +# per-column water bottom (fine-grid index), below-WB mask, frozen water field, and a linear v(z) trend +wbcol = [clamp(round(Int, wbdepth(ix)), 1, nz) for ix = 1:nx] +vtrue = true_model() +below = [iz >= wbcol[ix] for iz = 1:nz, ix = 1:nx] +water_field = vwater .* .!below +a0, b0 = linear_vtrend(vtrue, wbcol) +vtrend = iz -> a0 + b0 * (iz - 1) + +# build the RBF parameterization at a given frequency: tapered node cloud -> operator -> fit coefficients. +# precondition = false keeps the raw partition-of-unity kernel (best for a least-squares FIT / representation +# demo; the column-balancing precondition=true default is for reduced-parameterization inversion instead). +function build_rbf(freq) + nodes, deltas = tapered_rbf_nodes(nz, nx; wb = wbcol, vtrend = vtrend, freq = freq, dz = dz, + ppw_top = ppw_top, ppw_bot = ppw_bot, support = support) + M = size(nodes, 2) + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, nz, nx), nodes; delta = deltas, precondition = false) + P = JopDiagonal(Float64.(below)) ∘ A # freeze the water column + pou = A * ones(domain(A)) + c = convert(Matrix, P) \ vec(vtrue .* below) # fit coefficients to the true sediment + vmodel = reshape(P * c, nz, nx) .+ water_field + relerr = (vmodel .- vtrue) ./ vtrue + rms = sqrt(sum(abs2, relerr[below]) / count(below)) + @info "$(flabel(freq)) Hz: $M nodes ($(round(Int, count(below) / M))x reduction), min coverage $(round(minimum(pou[below]), digits=3)), RMS $(round(rms, digits=4))" + (; nodes, vmodel, relerr, rms, pou, M) +end + +res = [build_rbf(f) for f in freqs] + +# --------------------------------------------------------------------------- +# figure: 4 rows x 2 cols +# true model | coverage ; 6 Hz RBF | 6 Hz error ; 3 Hz | error ; 1.5 Hz | error +# --------------------------------------------------------------------------- +ext = [0.5, nx + 0.5, nz + 0.5, 0.5] +wbline = [wbdepth(ix) for ix = 1:nx] +vmin, vmax = 1500.0, maximum(vtrue) +emax = maximum(maximum(abs, r.relerr[below]) for r in res) # common error color scale + +function modelpanel(pos, img, ttl; nodes = nothing) + subplot(4, 2, pos) + im = imshow(img, extent = ext, aspect = "auto", cmap = "viridis", vmin = vmin, vmax = vmax) + plot(1:nx, wbline, "w-", lw = 1.2) + nodes === nothing || scatter(nodes[2, :], nodes[1, :], s = node_ms, c = "white", edgecolors = "k", linewidths = 0.3) + xlim(1, nx); ylim(nz, 1); xlabel("x"); ylabel("z"); title(ttl); colorbar(im, fraction = 0.046) +end + +function errpanel(pos, relerr, ttl) + subplot(4, 2, pos) + img = copy(relerr); img[.!below] .= NaN + im = imshow(img, extent = ext, aspect = "auto", cmap = "seismic", vmin = -emax, vmax = emax) + plot(1:nx, wbline, "k-", lw = 0.8) + xlim(1, nx); ylim(nz, 1); xlabel("x"); ylabel("z"); title(ttl) + cb = colorbar(im, fraction = 0.046); cb.set_label("Δv / v") +end + +figure(figsize = (12, 16)) + +# row 1: true model | coverage (finest frequency) +modelpanel(1, vtrue, "true model") +subplot(4, 2, 2) +covimg = copy(res[1].pou); covimg[.!below] .= NaN +imc = imshow(covimg, extent = ext, aspect = "auto", cmap = "magma", vmin = 0.999, vmax = 1.001) +plot(1:nx, wbline, "c-", lw = 1.2); xlim(1, nx); ylim(nz, 1); xlabel("x"); ylabel("z") +title("coverage below WB A·1 ≈ 1 ($(flabel(freqs[1])) Hz, water frozen white)"); colorbar(imc, fraction = 0.046) + +# rows 2-4: one frequency per row +for (i, f) in enumerate(freqs) + r = res[i] + modelpanel(2i + 1, r.vmodel, "$(flabel(f)) Hz RBF ($(r.M) nodes, ppw $(flabel(ppw_top))→$(flabel(ppw_bot)))"; nodes = r.nodes) + errpanel(2i + 2, r.relerr, @sprintf("%s Hz error (RMS %.4f)", flabel(f), r.rms)) +end + +tight_layout() +savefig(joinpath(OUT, "jopRBF_waterbottom_2d.png"), dpi = 120); close() +println("wrote ", joinpath(OUT, "jopRBF_waterbottom_2d.png")) diff --git a/docs/JopRBF/Project.toml b/docs/JopRBF/Project.toml new file mode 100644 index 0000000..88538ca --- /dev/null +++ b/docs/JopRBF/Project.toml @@ -0,0 +1,7 @@ +[deps] +JetPack = "24ef3835-3876-54c3-8a7a-956cf69ca0b2" +Jets = "2a57b368-ab28-5ba9-84aa-637fe7991822" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/JopRBF/images/jopRBF_basis_2d.png b/docs/JopRBF/images/jopRBF_basis_2d.png new file mode 100644 index 0000000..2463bd4 Binary files /dev/null and b/docs/JopRBF/images/jopRBF_basis_2d.png differ diff --git a/docs/JopRBF/images/jopRBF_smoothness_2d.png b/docs/JopRBF/images/jopRBF_smoothness_2d.png new file mode 100644 index 0000000..3f6899d Binary files /dev/null and b/docs/JopRBF/images/jopRBF_smoothness_2d.png differ diff --git a/docs/JopRBF/images/jopRBF_waterbottom_2d.png b/docs/JopRBF/images/jopRBF_waterbottom_2d.png new file mode 100644 index 0000000..babd6f6 Binary files /dev/null and b/docs/JopRBF/images/jopRBF_waterbottom_2d.png differ diff --git a/docs/JopRBF/rbf_tapered_nodes.jl b/docs/JopRBF/rbf_tapered_nodes.jl new file mode 100644 index 0000000..70762bf --- /dev/null +++ b/docs/JopRBF/rbf_tapered_nodes.jl @@ -0,0 +1,104 @@ +#= +rbf_tapered_nodes.jl -- generate a depth-tapered scattered node cloud for JopRBF, with NO external mesher. +Pure base Julia. Rows CONFORM to the water bottom: each row sits at a fixed offset from the (per-column) WB +and the offsets march downward by the local spacing, so the cloud is dense near the water bottom and coarser +at depth (no horizontal rows cross-cutting a dipping WB). Node spacing follows a points-per-wavelength rule, +spacing(z) = vtrend(z) / (freq * ppw(z)), with ppw LINEAR in depth from `ppw_top` at the water bottom to +`ppw_bot` at the model bottom. Rows are brick-offset (every other row shifted half a spacing) for near- +hexagonal packing; a few rows sit just above the WB (a capped "ghost" collar), a flat collar sits just below +the model bottom, and the lateral edges x=1/x=nx are pinned in every row, so all four boundaries get two- +sided support and the fit does not overshoot at any edge. +=# + +""" + nodes, deltas = tapered_rbf_nodes(nz, nx; wb, vtrend, freq, dz, ppw_top=5.0, ppw_bot=2.0, support=1.6, nghost=1, ghost_cap=4.0, nghost_bot=1) + +Depth-tapered scattered RBF node cloud on an `nz x nx` grid, in fine-grid index coordinates, with rows that +CONFORM to (run parallel to) the water bottom. `wb` is a per-column `Vector` of fine-grid depth indices +(length `nx`) or a scalar; `vtrend(iz)` returns a reference velocity [m/s] at depth index `iz`; `dz` is the +grid spacing [m]; `support` sets each node's compact support radius as a multiple of its local spacing. Each +row sits at a fixed offset from the (per-column) water bottom (`z = wb[ix] + offset`); the offsets march +downward by the local spacing `vtrend(z) / (freq * ppw(z)) / dz`, so the cloud is dense at the water bottom +and coarser at depth (ppw LINEAR in depth from `ppw_top` at the WB to `ppw_bot` at the model bottom). +`nghost` rows sit just ABOVE the water bottom (each offset capped at `ghost_cap` fine samples, so they stay +near the WB even at low frequency) to give the near-WB fit two-sided support and prevent boundary overshoot. +`nghost_bot` FLAT rows sit just below the model bottom (mirror of the WB ghost) and the lateral edges +x=1/x=nx are pinned in every row, so all four boundaries get two-sided support. Rows are brick-offset (every +other row shifted half a spacing) for near-hexagonal packing, clipped to the grid. Returns `nodes` (2 x M, +row 1 = z, row 2 = x) and `deltas` (2 x M per-node isotropic support radii), both in fine-grid index units, +ready for `JopRBF(...; delta=deltas)`. +""" +function tapered_rbf_nodes(nz::Integer, nx::Integer; wb, vtrend, freq::Real, dz::Real, + ppw_top::Real = 5.0, ppw_bot::Real = 2.0, support::Real = 1.6, + nghost::Integer = 1, ghost_cap::Real = 4.0, nghost_bot::Integer = 1) + wbcol = wb isa AbstractVector ? collect(Int, wb) : fill(Int(wb), nx) + length(wbcol) == nx || error("tapered_rbf_nodes -- wb vector must have length nx = $nx") + wbmin = minimum(wbcol) + ppw(z) = ppw_top + (ppw_bot - ppw_top) * clamp((z - (wbmin - 1)) / (nz - wbmin), 0.0, 1.0) + spacing(z) = vtrend(z) / (freq * ppw(z)) / dz # node spacing in FINE SAMPLES at depth index z + + # row offsets from the water bottom (fine samples): `nghost` capped rows ABOVE the WB (two-sided support), + # the WB row (offset 0), then march DOWN by the local spacing to the model bottom. + offs = Float64[] + o = 0.0 + for _ = 1:nghost + s = spacing(max(wbmin + o, 1.0)); s > 0 || error("tapered_rbf_nodes -- spacing must be positive; got $s") + o -= min(s, ghost_cap); pushfirst!(offs, o) + end + push!(offs, 0.0) + o = 0.0 + while wbmin + o <= nz + s = spacing(wbmin + o); s > 0 || error("tapered_rbf_nodes -- spacing must be positive; got $s") + o += s; wbmin + o <= nz && push!(offs, o) + end + + # brick-offset lateral positions for a row of spacing s, always pinning the x=1 and x=nx boundaries so the + # left/right edges get two-sided support (no one-sided lateral-edge overshoot). + rowx(s, k) = begin + xs = Float64[]; x = isodd(k) ? 1.0 : 1.0 + 0.5s + while x <= nx; push!(xs, x); x += s; end + (isempty(xs) || xs[1] > 1.0) && pushfirst!(xs, 1.0) + xs[end] < nx && push!(xs, Float64(nx)); xs + end + + # place brick-offset rows PARALLEL to the water bottom: node depth = wbcol[ix] + offset, clipped to grid. + zp = Float64[]; xp = Float64[]; rp = Float64[] + for (k, off) in enumerate(offs) + s = spacing(wbmin + max(off, 0.0)) # lateral spacing at this row's reference depth + for x in rowx(s, k) + ix = clamp(round(Int, x), 1, nx); z = wbcol[ix] + off + 1.0 <= z <= nz && (push!(zp, z); push!(xp, x); push!(rp, support * s)) + end + end + + # bottom collar: `nghost_bot` FLAT rows just BELOW the model bottom (mirror of the WB ghost) so the deepest + # cells get two-sided support instead of a one-sided edge (the model bottom is flat, so these rows are flat). + s_bot = spacing(nz); zb = float(nz) + for g = 1:nghost_bot + zb += min(s_bot, ghost_cap) + for x in rowx(s_bot, g + length(offs)) + push!(zp, zb); push!(xp, x); push!(rp, support * s_bot) + end + end + + nodes = permutedims(hcat(zp, xp)) # (2, M): row 1 = z, row 2 = x + nodes, permutedims(hcat(rp, rp)) # per-node isotropic support radius (2, M) +end + +""" + a, b = linear_vtrend(vmodel, wbcol) + +Least-squares linear fit `v(iz) = a + b*(iz-1)` [m/s] of the lateral-mean velocity below the water bottom, +a convenient `vtrend` for `tapered_rbf_nodes`. `wbcol` is the per-column water-bottom index. +""" +function linear_vtrend(vmodel::AbstractMatrix, wbcol::AbstractVector) + nz, nx = size(vmodel) + zc = Float64[]; vc = Float64[] + for iz = 1:nz + idx = [ix for ix = 1:nx if iz >= wbcol[ix]] + isempty(idx) && continue + push!(zc, iz - 1); push!(vc, sum(@view vmodel[iz, idx]) / length(idx)) + end + a, b = hcat(ones(length(zc)), zc) \ vc + a, b +end diff --git a/src/JetPack.jl b/src/JetPack.jl index 9987830..0d4056b 100644 --- a/src/JetPack.jl +++ b/src/JetPack.jl @@ -44,6 +44,7 @@ include("jop_taper.jl") include("jop_tanh.jl") include("jop_translation.jl") include("jop_CubicSpline.jl") +include("jop_RBF.jl") ###### JetPack with Flux extension diff --git a/src/jop_RBF.jl b/src/jop_RBF.jl new file mode 100644 index 0000000..ee5f81c --- /dev/null +++ b/src/jop_RBF.jl @@ -0,0 +1,277 @@ +""" + A = JopRBF(dom, rng, nodes; delta) + +Meshless interpolation from a *scattered* set of control nodes to a regular fine +grid, using **normalized compact-support radial basis functions** (Wendland C2). + +This is a reduced (coarse) model parameterization: the domain coefficients live on +an arbitrary cloud of points (not a grid), and the operator evaluates a smooth +field from them on the range grid. It is built for very large fine grids (e.g. +`1000^3`): every fine point depends on only the O(1) nodes within a support +radius, so the forward and adjoint are applied matrix-free with a node bucket +index; nothing dense is formed or stored. + +The field is the normalized (Shepard) combination +``` + Σ_j φ(r_j(x)) c_j +d(x) = ─────────────────────, r_j(x) = sqrt( Σ_k ((x_k - ξ_{j,k}) / δ_{k,j})^2 ) + Σ_j φ(r_j(x)) +``` +with the Wendland C2 kernel `φ(r) = (1-r)^4 (4r+1)` on `0 ≤ r < 1` (zero beyond). +Properties, all exercised by the unit tests: + +* **Compact support / scalable** -- `φ` vanishes for `r ≥ 1`, so node `j` only + influences fine points inside the ellipsoid of radii `δ_{:,j}` about `ξ_j`. +* **C2 smooth** -- `φ`, `φ'` and `φ''` vanish at `r = 1`, so the field is twice + continuously differentiable (no bilinear-style creases / hot zones). +* **Reproduces constants** -- the normalization is a partition of unity, so a + constant coefficient vector maps to that constant wherever there is coverage. +* **No overshoot** -- `φ ≥ 0`, so the field is a convex combination of nearby + coefficients and stays within their min/max (no ringing / halos). + +# Arguments + +`dom`: operator domain, `JetSpace(T, M)`, the `M` scattered node coefficients.\\ +`rng`: operator range, `JetSpace(T, n1[, n2[, n3]])`, the fine grid. The number of +active spatial dimensions is `ndims(rng)` (1, 2 or 3).\\ +`nodes`: a `ndims(rng) x M` matrix of node coordinates in *fine-grid index +coordinates* (fine sample `i` along an axis sits at coordinate `i`, so +coordinates lie in `[1, n_k]`). Placing some nodes on the boundary of the model +guarantees every fine point is covered (no extrapolation holes).\\ +`delta`: the compact-support radius, in fine-grid index units. It may be + +* a **scalar** (isotropic, the same radius for every node and axis); +* a length-`ndims(rng)` **vector** (per axis, e.g. for anisotropic grid spacing); + or +* a `ndims(rng) x M` **matrix** (per node and axis). This is the multiresolution + form: each node `j` carries its own support radii `δ_{:,j}`, so a + points-per-wavelength node cloud (dense with small support near the surface, + sparse with large support at depth) resolves fine detail where nodes are dense + and stays smooth where they are sparse. + +A fine point is influenced by node `j` when it is within `δ_{:,j}` of `ξ_j` along +every axis. The support must be large enough that every fine point of interest is +within some node's support; fine points not covered by any node (`Σ_j φ = 0`) map +to zero. + +`precondition` (default `true`): scale each node coefficient by `1/‖A e_j‖` (the +column norm, computed exactly from the kernel at construction) so the operator's +columns are ~equal norm and a unit coefficient maps to an order-unity model +perturbation. This removes the node support-size bias (dense small-support nodes vs +sparse large-support nodes), which is what a reduced-parameterization inverse +problem wants: a diagonal (Jacobi) preconditioner that makes the reduced +Gauss-Newton Gram `AᵀA` have ~unit diagonal. It trades away the raw kernel's exact +partition-of-unity / constant-reproduction (the all-ones coefficient no longer maps +to a constant); set `precondition = false` for the pure normalized RBF. +`precondition_weight` (optional, a range-sized array `w`): balance `‖diag(w) A e_j‖` +instead of `‖A e_j‖`, for when a fine-grid diagonal `diag(w)` is applied downstream +(e.g. an FWI illumination), so the FULL operator `diag(w) A` is column-balanced. + +# Examples + +```julia +# 2D: 300 scattered nodes onto a 400 x 500 grid +nz, nx, M = 400, 500, 300 +nodes = vcat((1 .+ (nz-1).*rand(1,M)), (1 .+ (nx-1).*rand(1,M))) # 2 x M, index coords +A = JopRBF(JetSpace(Float64,M), JetSpace(Float64,nz,nx), nodes; delta = 60.0) +c = rand(domain(A)) +d = A*c +``` + +```julia +# per-node (multiresolution) support: each node gets its own radii (2 x M) +deltas = vcat(fill(15.0, 1, M), fill(15.0, 1, M)) # here isotropic, but may vary per node +A = JopRBF(JetSpace(Float64,M), JetSpace(Float64,nz,nx), nodes; delta = deltas) +``` +""" +function JopRBF(dom::JetSpace{T,1}, rng::JetSpace{T,D}, nodes::AbstractMatrix; delta, + precondition::Bool = true, precondition_weight = nothing) where {T,D} + (D >= 1 && D <= 3) || error("JopRBF supports 1, 2 or 3 spatial (range) dimensions, got $D") + M = size(dom)[1] + n = size(rng) + size(nodes) == (D, M) || + error("nodes must be (ndims(range)=$D) x (length(domain)=$M), got $(size(nodes))") + precondition_weight === nothing || size(precondition_weight) == n || + error("precondition_weight must match the range size $(n), got $(size(precondition_weight))") + + δ = _delta_matrix(delta, D, M) # (D, M) per-node per-axis radii + all(x -> x > 0, δ) || error("delta must be positive") + + # bucket size = the largest support radius per axis, so that any node covering + # a fine point sits within one bucket of it (exact neighbor search). + bucketsize = ntuple(k -> maximum(@view δ[k, :]), D) + + nodesf = Matrix{Float64}(nodes) + buckets = _build_buckets(nodesf, bucketsize) + invw = _build_invw(n, nodesf, δ, bucketsize, buckets) + # optional construction-time diagonal (Jacobi) preconditioner: scale each node coefficient by 1/‖A e_j‖ + # so the operator's columns are ~equal norm (a unit coefficient -> an order-unity model perturbation, + # removing the node support-size bias). Computed EXACTLY from the kernel; ones when precondition=false. + colscale = precondition ? _column_scale(n, nodesf, δ, invw, precondition_weight) : ones(Float64, M) + + JopLn(; dom, rng, df! = JopRBF_df!, df′! = JopRBF_df′!, + s = (; nodes = nodesf, delta = δ, bucketsize, buckets, invw, colscale)) +end + +export JopRBF + +_delta_matrix(delta::Number, D, M) = fill(Float64(delta), D, M) +function _delta_matrix(delta::AbstractVector, D, M) + length(delta) == D || error("delta vector length ($(length(delta))) must equal number of range dimensions ($D)") + Float64[delta[k] for k = 1:D, j = 1:M] +end +function _delta_matrix(delta::AbstractMatrix, D, M) + size(delta) == (D, M) || error("delta matrix must be (ndims(range)=$D) x (length(domain)=$M), got $(size(delta))") + Matrix{Float64}(delta) +end + +""" +Wendland C2 radial basis function: `φ(r) = (1-r)^4 (4r+1)` for `0 ≤ r < 1`, else +`0`. Positive definite in up to three dimensions, non-negative, and C2 (value, +first and second derivatives vanish at `r = 1`). +""" +@inline function wendland_c2(r::Float64) + if r < 1.0 + om = 1.0 - r + return om * om * om * om * (4.0 * r + 1.0) + else + return 0.0 + end +end + +# assign each node to a bucket of size `bucketsize` (per axis) for O(1) neighbor search +function _build_buckets(nodes::Matrix{Float64}, bucketsize::NTuple{D,Float64}) where {D} + M = size(nodes, 2) + buckets = Dict{NTuple{D,Int}, Vector{Int}}() + for j = 1:M + b = ntuple(k -> floor(Int, nodes[k, j] / bucketsize[k]), D) + push!(get!(() -> Int[], buckets, b), j) + end + buckets +end + +# precompute the normalization: invw[x] = 1 / Σ_j φ(r_j(x)) (0 where uncovered). +# Built once from node positions and reused for every apply. +function _build_invw(n::NTuple{D,Int}, nodes::Matrix{Float64}, delta::Matrix{Float64}, bucketsize::NTuple{D,Float64}, buckets) where {D} + invw = zeros(Float64, n) + Threads.@threads :static for I in CartesianIndices(invw) + x = Tuple(I) + w = _accumulate_weight(x, nodes, delta, bucketsize, buckets) + @inbounds invw[I] = w > 0.0 ? 1.0 / w : 0.0 + end + invw +end + +# sum of kernel weights at fine coordinate `x` over nodes in the 3^D neighbor buckets +@inline function _accumulate_weight(x::NTuple{D,Int}, nodes::Matrix{Float64}, delta::Matrix{Float64}, bucketsize::NTuple{D,Float64}, buckets) where {D} + bx = ntuple(k -> floor(Int, x[k] / bucketsize[k]), D) + w = 0.0 + for off in CartesianIndices(ntuple(_ -> -1:1, D)) + ob = Tuple(off) + b = ntuple(k -> bx[k] + ob[k], D) + nb = get(buckets, b, nothing) + nb === nothing && continue + for j in nb + r2 = 0.0 + @inbounds for k = 1:D + dk = (x[k] - nodes[k, j]) / delta[k, j] + r2 += dk * dk + end + if r2 < 1.0 + w += wendland_c2(sqrt(r2)) + end + end + end + w +end + +# exact per-node column norm of the operator (optionally weighted): for node j, sum over its support box the +# squared normalized-basis value (φ_j·invw) times the optional fine-grid weight². Enumerates the SAME (node, +# fine-point) pairs as the adjoint, so it is exact. Returns the preconditioner scale = median(‖·‖)/‖A e_j‖, +# clamped to [0.1, 10], so a unit coefficient maps to an order-unity model perturbation. +_median(v::AbstractVector{Float64}) = (s = sort(v); m = length(s); iseven(m) ? (s[m ÷ 2] + s[m ÷ 2 + 1]) / 2 : s[(m + 1) ÷ 2]) + +function _column_scale(n::NTuple{D,Int}, nodes::Matrix{Float64}, delta::Matrix{Float64}, invw, weight) where {D} + M = size(nodes, 2) + d = zeros(Float64, M) + Threads.@threads :static for j = 1:M + xlo = ntuple(k -> max(1, ceil(Int, nodes[k, j] - delta[k, j])), D) + xhi = ntuple(k -> min(n[k], floor(Int, nodes[k, j] + delta[k, j])), D) + acc = 0.0 + for I in CartesianIndices(ntuple(k -> xlo[k]:xhi[k], D)) + x = Tuple(I) + r2 = 0.0 + @inbounds for k = 1:D + dk = (x[k] - nodes[k, j]) / delta[k, j] + r2 += dk * dk + end + if r2 < 1.0 + w = weight === nothing ? 1.0 : Float64(weight[I]) + a = wendland_c2(sqrt(r2)) * invw[I] * w + acc += a * a + end + end + @inbounds d[j] = sqrt(acc) + end + pos = filter(>(0), d) + dref = isempty(pos) ? 1.0 : _median(pos) + [clamp(dref / max(dj, dref * 1e-3), 0.1, 10.0) for dj in d] +end + +# +# forward: thread over fine points, gather nearby nodes via the bucket index +# (each fine point writes only its own output, so no data races). +# +function JopRBF_df!(d::AbstractArray{T,D}, c::AbstractVector{T}; nodes, delta, bucketsize, buckets, invw, colscale, kwargs...) where {T,D} + cc = colscale .* c # construction-time column preconditioner (ones if off) + Threads.@threads :static for I in CartesianIndices(d) + x = Tuple(I) + bx = ntuple(k -> floor(Int, x[k] / bucketsize[k]), D) + num = 0.0 + for off in CartesianIndices(ntuple(_ -> -1:1, D)) + ob = Tuple(off) + b = ntuple(k -> bx[k] + ob[k], D) + nb = get(buckets, b, nothing) + nb === nothing && continue + for j in nb + r2 = 0.0 + @inbounds for k = 1:D + dk = (x[k] - nodes[k, j]) / delta[k, j] + r2 += dk * dk + end + if r2 < 1.0 + @inbounds num += wendland_c2(sqrt(r2)) * cc[j] + end + end + end + @inbounds d[I] = convert(T, num * invw[I]) + end + d +end + +# +# adjoint: thread over nodes, iterate each node's fine-grid bounding box (each +# node writes only its own coefficient, so no data races). This enumerates the +# same (fine point, node) pairs as the forward, so it is the exact transpose. +# +function JopRBF_df′!(c::AbstractVector{T}, d::AbstractArray{T,D}; nodes, delta, bucketsize, buckets, invw, colscale, kwargs...) where {T,D} + n = size(d) + Threads.@threads :static for j = 1:length(c) + xlo = ntuple(k -> max(1, ceil(Int, nodes[k, j] - delta[k, j])), D) + xhi = ntuple(k -> min(n[k], floor(Int, nodes[k, j] + delta[k, j])), D) + acc = 0.0 + for I in CartesianIndices(ntuple(k -> xlo[k]:xhi[k], D)) + x = Tuple(I) + r2 = 0.0 + @inbounds for k = 1:D + dk = (x[k] - nodes[k, j]) / delta[k, j] + r2 += dk * dk + end + if r2 < 1.0 + @inbounds acc += wendland_c2(sqrt(r2)) * (convert(Float64, d[I]) * invw[I]) + end + end + @inbounds c[j] = convert(T, acc * colscale[j]) + end + c +end diff --git a/test/jop_RBF.jl b/test/jop_RBF.jl new file mode 100644 index 0000000..574fd48 --- /dev/null +++ b/test/jop_RBF.jl @@ -0,0 +1,166 @@ +using Jets, LinearAlgebra, Test, JetPack + +# a regular lattice of nodes (in fine-index coordinates) that includes the domain +# boundary, optionally jittered, guaranteeing full coverage for delta >= 1.5*h. +function lattice_nodes(n::NTuple{D,Int}, h::Real; jitter = 0.0, rng = nothing) where {D} + axes1d = ntuple(k -> unique(vcat(collect(1.0:h:n[k]), Float64(n[k]))), D) + grid = Iterators.product(axes1d...) + pts = [collect(Float64, p) for p in grid] + nodes = reduce(hcat, pts) + if jitter > 0 + r = rng === nothing ? () : (rng,) + for j = 1:size(nodes, 2), k = 1:D + interior = nodes[k, j] > 1.5 && nodes[k, j] < n[k] - 0.5 + interior && (nodes[k, j] += jitter * (rand(r...) - 0.5)) + end + end + nodes +end + +@testset "JopRBF, wendland_c2 kernel is C2 and non-negative" begin + φ = JetPack.wendland_c2 + @test φ(0.0) == 1.0 + @test φ(1.0) == 0.0 + @test φ(1.5) == 0.0 + @test all(φ(r) >= 0 for r in range(0, 1, length = 101)) + # value, first and second derivative all vanish at r -> 1^- (=> C2 at support edge) + h = 1e-4 + r0 = 1.0 - h + d1 = (φ(r0 + h/2) - φ(r0 - h/2)) / h + d2 = (φ(r0 + h) - 2φ(r0) + φ(r0 - h)) / h^2 + @test isapprox(φ(r0), 0.0; atol = 1e-8) + @test isapprox(d1, 0.0; atol = 1e-3) + @test isapprox(d2, 0.0; atol = 1e-1) +end + +@testset "JopRBF, dot product test - $(D)D" for D in (1, 2, 3) + n = D == 1 ? (151,) : D == 2 ? (61, 47) : (23, 19, 17) + M = 40 + nodes = zeros(D, M) + for k = 1:D, j = 1:M + nodes[k, j] = 1 + (n[k] - 1) * rand() + end + delta = 0.5 * minimum(n) + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = delta) + lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A))) + @test lhs ≈ rhs +end + +@testset "JopRBF, reproduces constants (partition of unity) - $(D)D" for D in (1, 2, 3) + n = D == 1 ? (120,) : D == 2 ? (60, 50) : (24, 20, 18) + h = D == 1 ? 12 : D == 2 ? 12 : 6 + nodes = lattice_nodes(n, h) + M = size(nodes, 2) + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = 1.5h, precondition = false) + d = A * ones(domain(A)) + @test all(x -> isapprox(x, 1.0; atol = 1e-12), d) # full coverage => constant reproduced everywhere +end + +@testset "JopRBF, no overshoot (convex combination) - $(D)D" for D in (2, 3) + n = D == 2 ? (60, 50) : (24, 20, 18) + h = D == 2 ? 12 : 6 + nodes = lattice_nodes(n, h) + M = size(nodes, 2) + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = 1.5h, precondition = false) + c = rand(domain(A)) # coefficients in [0,1] + d = A * c + @test all(x -> -1e-12 <= x <= 1 + 1e-12, d) + @test minimum(d) >= minimum(c) - 1e-12 + @test maximum(d) <= maximum(c) + 1e-12 +end + +@testset "JopRBF, compact support - influence vanishes beyond delta" begin + n = (80, 90) + nodes = [20.0 60.0; # node 1 at (20,30), node 2 at (60,70) + 30.0 70.0] + delta = 18.0 + A = JopRBF(JetSpace(Float64, 2), JetSpace(Float64, n...), nodes; delta = delta) + e1 = zeros(domain(A)); e1[1] = 1.0 + d = A * e1 + for I in CartesianIndices(d) + x = Tuple(I) + r = sqrt((x[1] - nodes[1,1])^2 + (x[2] - nodes[2,1])^2) / delta + if r >= 1.0 + @test d[I] == 0.0 # zero outside node 1's support + end + end + @test any(d .> 0) # nonzero somewhere inside the support +end + +@testset "JopRBF, anisotropic (per-axis) delta - dot product test" begin + n = (50, 40) + M = 30 + nodes = zeros(2, M) + for j = 1:M + nodes[1, j] = 1 + (n[1] - 1) * rand() + nodes[2, j] = 1 + (n[2] - 1) * rand() + end + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = [10.0, 20.0]) + lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A))) + @test lhs ≈ rhs +end + +@testset "JopRBF, per-node (multiresolution) delta" begin + # each node carries its own support radii (2 x M). Small support near the + # surface, large support at depth (a points-per-wavelength style cloud). + n = (120, 90) + nodes = lattice_nodes(n, 12) + M = size(nodes, 2) + # per-node radius grows with depth (row-1 coordinate) + deltas = zeros(2, M) + for j = 1:M + z = nodes[1, j] + rad = 14.0 + 40.0 * (z - 1) / (n[1] - 1) # 14 shallow -> 54 deep + deltas[1, j] = rad + deltas[2, j] = rad + end + A = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = deltas, precondition = false) + lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A))) + @test lhs ≈ rhs + # constants still reproduced where covered + d = A * ones(domain(A)) + covered = d .> 0 + @test all(x -> isapprox(x, 1.0; atol = 1e-12), d[covered]) + @test count(covered) == length(d) # full coverage for these radii + # a shallow node influences a smaller region than a deep node + jshallow = argmin(nodes[1, :]); jdeep = argmax(nodes[1, :]) + supp(j) = (e = zeros(domain(A)); e[j] = 1.0; count(!iszero, A * e)) + @test supp(jdeep) > supp(jshallow) +end + +@testset "JopRBF, single precision support" begin + n = (60, 50) + nodes = lattice_nodes(n, 12) + M = size(nodes, 2) + A = JopRBF(JetSpace(Float32, M), JetSpace(Float32, n...), Float32.(nodes); delta = 18.0f0, precondition = false) + @test eltype(domain(A)) == Float32 + lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A))) + @test lhs ≈ rhs + d = A * ones(domain(A)) + @test all(x -> isapprox(x, 1f0; atol = 1f-5), d) +end + +@testset "JopRBF, precondition=true equalizes column norms (multiresolution)" begin + # a node cloud with support that varies ~4x with depth => the raw kernel has a big column-norm (support- + # size) imbalance; the default precondition=true rescales each coefficient by 1/‖A e_j‖ so the columns are + # ~equal norm (order unity), while keeping an exact adjoint. This is the Jacobi preconditioner for a + # reduced-parameterization inverse problem, moved into the operator at construction. + n = (120, 90) + nodes = lattice_nodes(n, 12) + M = size(nodes, 2) + deltas = zeros(2, M) + for j = 1:M + rad = 14.0 + 40.0 * (nodes[1, j] - 1) / (n[1] - 1) # 14 shallow -> 54 deep + deltas[1, j] = rad; deltas[2, j] = rad + end + colnorms(A) = [norm(A * setindex!(zeros(domain(A)), 1.0, j)) for j in 1:M] + spread(v) = maximum(v) / minimum(v) + A0 = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = deltas, precondition = false) + A1 = JopRBF(JetSpace(Float64, M), JetSpace(Float64, n...), nodes; delta = deltas, precondition = true) + @test spread(colnorms(A0)) > 3.0 # raw kernel: large support-size imbalance + @test spread(colnorms(A1)) < 1.5 # preconditioned: ~order unity + lhs, rhs = dot_product_test(A1, rand(domain(A1)), rand(range(A1))) + @test lhs ≈ rhs # exact adjoint preserved under the preconditioner +end + +nothing diff --git a/test/runtests.jl b/test/runtests.jl index f7fd9ae..5c78d00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,7 @@ for filename in ( "jop_taper.jl", "jop_tanh.jl", "jop_translation.jl", - "jop_CubicSpline.jl") + "jop_CubicSpline.jl", + "jop_RBF.jl") include(filename) end