Add JopRBF: scattered-node compact-support RBF model parameterization - #48
Add JopRBF: scattered-node compact-support RBF model parameterization#48jkwashbourne-oss wants to merge 6 commits into
JopRBF: scattered-node compact-support RBF model parameterization#48Conversation
- JopRBF operator: normalized compact-support Wendland-C2 RBF from a scattered node cloud to a regular fine grid; scalar/per-axis/per-node support radii; matrix-free forward+adjoint via a node bucket index (scales to large grids). - Unit tests: kernel C2, dot-product (adjoint) 1D/2D/3D, partition of unity, compact support, no overshoot, anisotropic and per-node delta, Float32. - docs/JopRBF: water-bottom freeze-workflow demo (Gmsh-meshed nodes, per-frequency comparison) + kernel/smoothness demo, explanatory JopRBF-demo.md, and a local project (Project.toml) so PyPlot/Gmsh stay out of JetPack.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #48 +/- ##
==========================================
+ Coverage 95.94% 96.04% +0.10%
==========================================
Files 40 41 +1
Lines 1086 1215 +129
==========================================
+ Hits 1042 1167 +125
- Misses 44 48 +4 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…operator at construction Adds precondition::Bool=true + optional precondition_weight to JopRBF. At construction it computes the exact per-node column norm (reusing the adjoint support-box enumeration) and folds colscale_j = clamp(median/norm, 0.1, 10) into the forward and adjoint, so the operator columns are ~equal norm: a unit coefficient maps to an order-unity model perturbation (a Jacobi preconditioner for a reduced-parameterization inverse problem; removes the node support-size bias), adjoint still exact. NOTE: precondition=true trades away the raw kernel exact partition-of-unity and constant-reproduction (the all-ones coefficient no longer maps to a constant); pass precondition=false for the pure kernel. Tests: partition-of-unity/constant/no-overshoot tests use precondition=false; added a precondition column-norm-equalization test; dot-product/compact-support tests validate the preconditioned adjoint at the default.
…cloud The water-bottom demo now places RBF centers with a pure base-Julia depth-tapered structured cloud (rbf_tapered_nodes.jl: brick-offset rows, spacing = v(z)/(freq*ppw(z)), ppw tapering with depth) instead of a Gmsh mesh. Removes Gmsh from docs/JopRBF/Project.toml and Manifest, rewrites JopRBF_waterbottom_demo.jl, and removes all Gmsh mentions from JopRBF-demo.md. Figure regenerated (6/3/1.5 Hz to 474/128/43 nodes).
…sided support tapered_rbf_nodes now lays nodes in rows PARALLEL to the water bottom (no horizontal rows cross-cutting a dipping WB), with a capped ghost row above the WB (nghost/ghost_cap), a flat collar below the model bottom (nghost_bot), and the lateral edges x=1/x=nx pinned in every row, so all four boundaries get two-sided support and the fit does not overshoot at any edge. Water-bottom demo: build the true model at the SAME integer WB as the below-WB mask (was fractional), which removes the spurious WB overshoot band (water velocity was leaking into the top sediment cell). Regenerated figure and refreshed the node-count/RMS numbers (6 Hz 608/0.0008, 3 Hz 186/0.0013, 1.5 Hz 65/0.0032, precondition=false so coverage reads ~1.0).
| 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 |
There was a problem hiding this comment.
precondition is missing on the call docstring on line 2
| 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‖` |
| 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 |
There was a problem hiding this comment.
is there any safeguard on delta? if it is too small, wouldn't that cause some output grid points to be "orphaned".
|
Very nice addition to JetPack! It will be interesting to compare the RBF with Bspline with the same number of degrees of freedom and same node locations. RBF certainly offers much more flexibility, but that comes at a cost; apart from being a heavier operator, my guess is that the RBF is not as well conditioned as the Bspline operator and leaves clear footprints of the node's locations. But of course, it is more powerful in enforcing sharp irregular boundaries and variable smoothing by regions. |
| 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 |
There was a problem hiding this comment.
this inversion \ could be costly in practice depending on which method Julia used to actually solve \. It would be nice to use a guaranteed matrix-free method such as lsqr or lsmr from IterativeSolvers or others, and run a reasonable number of iterations to see how fast and how well the RBF can be reconstructed. That also gives information about how well conditioned the RBF operator is.
|
I think the added docs would need to be tied into docs/src/index.html somehow. Up until now, we have only used doc strings. Perhaps the demo can go in the Examples repository instead? |
4 participants
Add
JopRBF: scattered-node compact-support RBF model parameterizationSummary
Adds a new linear operator,
JopRBF, that 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 forJetsinverse 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). No new dependencies are added to JetPack.The operator
δ. Forward and adjoint are matrix-free via a node bucket index (no denseN_fine x Mmatrix): 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 (fine point, node) pairs, so the adjoint is exact.φ,φ',φ''vanish atr = 1(no bilinear-style creases).φ ≥ 0, convex combination of nearby coefficients).δ_{k,j}matrix (multiresolution).1/Σφdepend only on node positions andδ; there is no linear solve.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
δ); no overshoot; anisotropic per-axis and per-nodeδ; Float32. Registered insrc/JetPack.jlandtest/runtests.jl. The full suite passes locally (Pkg.test()).Docs and demo (
docs/JopRBF/)JopRBF-demo.md: explains the operator, the water-bottom "freeze" workflow, and why normalized compact-support RBF is preferred over global RBF / triangulation-cubic / barycentric-linear.JopRBF_waterbottom_demo.jl: a 101x201 water-bottom example. The sediment below a sloping water bottom is meshed with Gmsh using a points-per-wavelength size field (spacing = λ(z)/ppw(z)), with the water bottom embedded so mesh nodes land on it; the water column is frozen viaP = S_below ∘ A. The figure compares the parameterization across inversion frequencies (lower frequency => longer wavelength => coarser mesh).JopRBF_demo.jl: kernel and C0-vs-C2 smoothness figures.docs/JopRBF/Project.toml, not in JetPack. Run withjulia --project=docs/JopRBF docs/JopRBF/JopRBF_waterbottom_demo.jl.Notes
JopRBFand its docs/tests.