\documentclass[11pt]{article} \usepackage[margin=1in]{geometry} \usepackage{amsmath,amssymb,amsthm,mathtools} \usepackage{bm} \usepackage{booktabs} \usepackage{graphicx} \usepackage{longtable} \usepackage{hyperref} \usepackage{enumitem} \usepackage{xcolor} \usepackage{float} \hypersetup{ colorlinks=true, linkcolor=blue!50!black, citecolor=blue!50!black, urlcolor=blue!50!black } \newtheorem{definition}{Definition} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newcommand{\R}{\mathbb{R}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\C}{\mathbb{C}} \newcommand{\dd}{\mathrm{d}} \newcommand{\calL}{\mathcal{L}} \newcommand{\calN}{\mathcal{N}} \newcommand{\bs}{\boldsymbol} \newcommand{\A}{\mathbf{A}} \newcommand{\G}{\mathbf{G}} \newcommand{\I}{\mathbf{I}} \newcommand{\cvec}{\mathbf{c}} \newcommand{\xvec}{\mathbf{x}} \newcommand{\yvec}{\mathbf{y}} \newcommand{\wvec}{\mathbf{w}} \newcommand{\zvec}{\mathbf{z}} \newcommand{\uvec}{\mathbf{u}} \DeclareMathOperator*{\argmin}{arg\,min} \title{Operator-Spline Neural Representations: \\ Continuous-Domain Operator Bases for Efficient Physical Learning} \author{OSNR Project Notes} \date{\today} \begin{document} \maketitle \begin{abstract} Operator-Spline Neural Representations study how known continuous structure can be represented and manipulated through discrete coefficients. For admissible constant-coefficient operators, polynomial and exponential B-splines provide compact local generators, operator-null-space reproduction and analog-to-digital filtering identities. Hermite generators provide explicit value and derivative coordinates. This document develops those continuous--discrete connections and records the implementation and empirical conditions under which they are useful. The central computational objects are continuous function, derivative and cross-basis Grams. Fixed generators and grids permit their reuse across coefficient updates; projection, differential energies and adjoints then become structured linear algebra. Periodic equal-spacing settings admit circulant or block-circulant realizations, whereas unequal-grid cross-Grams, finite-boundary operators and data-dependent statistics need their actual structure respected. The implementations address canonical endpoint handling, stable reciprocal-root factorizations, positive spectra, banded/bordered solves, exact polynomial-piece integration and local matrix/Horner execution. Grid refinement preserves a function only with the required subspace embedding; changing a feature family does not in general preserve sufficient statistics for unseen directions. Recent extensions apply this calculus to compact causal physical models. For fixed spline directions and an unclipped regression model, a stable first-order actuator permits a joint constrained quadratic fit after absorbing its pole-dependent scale into static spline coefficients. For monotone fitted laws, explicit saturation also admits finitely many ordered data regions with constrained quadratic subproblems; batched quadratic lower bounds prune most solves while retaining numerical optimality checks. This is not a claim that the clipped objective is globally convex or that every fit meets a fixed real-time bound. Combining bounded-function products with stable exponential response sums yields an exact infinite-response inner product and a compact source-bank descriptor. This is exponential response calculus around a cardinal cubic model, not an exponential B-spline neural activation experiment or a new universal-policy theorem. A thirty-pool temporal-task transfer experiment provides a further boundary: cardinal support planning solves 161/480 composed missions versus MLP 242/480, with sequential reuse matching both. Exact Bernstein specialization preserves audited development outputs while reducing planning time by $2.10$--$2.46\times$, but informative state acquisition and useful composition do not follow from correct local calculus. This experiment uses cubic fields, not the full exponential/Hermite physical-operator toolbox. A measured CT extension combines exact rectangular field queries, finite-domain hierarchical product matrices, incremental likelihood information and full Gaussian block-design solves. Its strongest classical control answers 47 of 72 late questions without confident mistakes, but fails the fixed coverage target and gains no measurement saving from targeted acquisition. Exact-output memory compilation does not by itself establish useful uncertainty or a spline-specific inspection advantage. The statistical and behavioral claims remain conditional. Generalized increments filter continuous innovations and need not be independent. Operator matching alone does not guarantee lower prediction risk or sparse support adaptation by ridge. Pooled fixed-feature statistics preserve a fitting objective, not every old prediction or private datum; protected supports and immutable archived programs provide different preservation contracts. Physical identification also does not imply useful closed-loop recovery: the accompanying controlled experiments include both gains and failed restoration or comparator criteria. The resulting foundation is an auditable set of representation, calculus and compilation mechanisms, with explicit boundaries on observability, model mismatch, numerical conditioning, memory and downstream performance. \end{abstract} \tableofcontents \section{Introduction} Coordinate networks have become a standard tool for representing continuous signals and fields. A typical implicit neural representation (INR) maps coordinates to field values through a multilayer perceptron, \[ f_\theta : \xvec \mapsto y, \] where the weights $\theta$ are optimized by stochastic gradient descent. SIRENs improve high-frequency representation by using sinusoidal activations, while PINNs add differential-equation residuals to the training objective. The OSNR thesis is that this approach is structurally misaligned for physical systems governed by known operators. It is inspired by the operator-based spline signal-processing program of Unser, Blu, Vetterli, and collaborators \cite{unser1993bspline1,unser1993bspline2,unser2005cardinal1,unser2005cardinal2,vetterli2002fri}. If a field is constrained by \[ L\{s\} = r, \] then the representation should be built from the operator $L$ itself. The role of learning or numerical inversion should be reduced to coefficient recovery in an already appropriate continuous function space. OSNR therefore replaces black-box nonlinear layers by operator-matched spline dictionaries. The spline coefficients play the role of the latent representation. The architecture is not a generic neural network with a different activation; it is a continuous-domain inverse-problem engine presented in neural-representation form, aligned with continuous-domain inverse-problem representer theorems \cite{gupta2018continuous,debarre2019hybrid,debarre2021composite}. The same principle also suggests a deterministic alternative to the trunk side of branch-trunk neural operator models such as DeepONet and MIONet \cite{lu2021deeponet,jin2022mionet}: the coordinate-to-field map can be a spline synthesis operator with exact calculus rather than a learned MLP. The empirical record tests this thesis rather than establishing universal superiority. \S\ref{sec:grown-topology} covers closed-form identification, grown-topology control and continual-learning comparisons; \S\ref{sec:matched-vs-sindy} compares operator-matched identification with SINDy and neural ODEs, including a negative blind-discovery control. \S\ref{sec:pde-vs-fno} studies PDE comparisons with neural operators, regime changes and partially known physics. Their measured gains depend on the information supplied, comparator, accuracy and computational accounting: knowing an operator family is not by itself a sufficient condition for an advantage. The later exact battery-response experiment, for example, admits a nearly equally fast training-free classical control; the real battery-aging fit loses to its matched polynomial comparator. \S\ref{sec:ssp-view} develops the sparse-stochastic-process connection under its stated assumptions, not a universal statistical dominance theorem. \S\ref{sec:rsi} records recursive-improvement experiments, while subsequent claim audits distinguish pooled fitting statistics, protected function supports, immutable archives and retained task performance. None alone establishes unrestricted self-improvement without forgetting or drift. The constructive result is a tested representation and calculus toolbox; independent-task benefit against strong controls remains the capability test. \section{Mathematical Background} \subsection{Cardinal spline representation} The classical cardinal spline model represents a continuous signal as \cite{unser1993bspline1,unser1993bspline2} \[ s(t) = \sum_{k \in \Z} c[k] \beta(t-k), \] where $\beta$ is a spline generator and $c[k]$ are discrete coefficients. The essential point is that a continuous function space is controlled by a discrete sequence. This makes exact digital processing of continuous objects possible, provided the generator is stable and the correct coefficient-domain filters are used. \subsection{Linear differential operators and null spaces} Let $L$ be a constant-coefficient differential operator with characteristic roots or poles \[ \bs{\alpha} = (\alpha_1,\ldots,\alpha_N). \] The null space is \[ \calN_L = \{s : L\{s\}=0\} = \mathrm{span}\{e^{\alpha_m t}\}_{m=1}^N \] with polynomial factors included for repeated roots. For the Helmholtz operator \[ L = \frac{\dd^2}{\dd x^2} + k^2, \] the poles are $\alpha=\pm j k$, and the null space consists of sinusoidal modes. \subsection{Cardinal exponential splines} Cardinal exponential splines are compactly supported spline generators matched to $\bs{\alpha}$ \cite{unser2005cardinal1}. Their Fourier-domain form is \[ \widehat{\beta}_{\bs{\alpha}}(\omega) = \prod_{m=1}^N \frac{1 - e^{\alpha_m - j\omega}}{j\omega-\alpha_m}. \] Integer shifts of $\beta_{\bs{\alpha}}$ generate a stable spline space under appropriate Riesz conditions. The basis reproduces exponential polynomials and supports exact operator calculus in the coefficient domain. \begin{definition}[Operator-matched OSNR dictionary] Given an operator $L$ with pole vector $\bs{\alpha}$ and physical knot spacing $T$, an OSNR dictionary is a matrix of samples \[ \A_{n,k} = \beta_{\bs{\alpha}}\!\left(\frac{x_n}{T}-k\right), \] augmented when necessary by explicit null-space columns. The continuous field is represented by \[ s(x_n) \approx (\A \cvec)_n. \] \end{definition} \subsection{Exponential Hermite splines} Cardinal E-splines attach the operator to a scalar coefficient stream. For higher-order boundary value problems, OSNR requires a vector-valued cardinal system whose coefficients store not only function values but also derivatives. The second-order exponential-polynomial Hermite generator of Schmitter, Badoual, Uhlmann, Fageot, and Unser \cite{schmitter2016hermite} provides exactly this structure. Let \[ \Phi(t)= \begin{bmatrix} \phi_0(t) & \phi_1(t) & \phi_2(t) \end{bmatrix}^{\top}. \] The Hermite spline expansion is \[ f(t) = \sum_{k\in\Z} \left( c_0[k]\phi_0(t-k) +c_1[k]\phi_1(t-k) +c_2[k]\phi_2(t-k) \right), \] with coefficient vectors \[ \mathbf{c}[k]= \begin{bmatrix} c_0[k] \\ c_1[k] \\ c_2[k] \end{bmatrix} = \begin{bmatrix} f(k) \\ f'(k) \\ f''(k) \end{bmatrix} \] for exactly interpolated data. The interpolation conditions are \[ \phi_p^{(r)}(k)=\delta_{pr}\delta_k, \qquad p,r\in\{0,1,2\},\quad k\in\Z. \] On the interval $[0,1]$, each channel has the exponential-polynomial form \[ \phi_p(t) = A_p+B_p t+C_p t^2+D_p t^3 +E_p e^{j\omega_0 t}+F_p e^{-j\omega_0 t}, \qquad \omega_0=\frac{2\pi}{M}. \] The constants $A_p,\ldots,F_p$ are determined by the six Hermite endpoint constraints \[ \phi_p^{(r)}(0)=\delta_{pr}, \qquad \phi_p^{(r)}(1)=0, \qquad r=0,1,2. \] The negative branch is fixed by the Hermite parity \[ \phi_0(-t)=\phi_0(t),\qquad \phi_1(-t)=-\phi_1(t),\qquad \phi_2(-t)=\phi_2(t), \] so the generator is compactly supported on $[-1,1]$. This construction is $C^2$, reproduces polynomials up to cubic degree, and reproduces the trigonometric modes $\sin(\omega_0 t)$ and $\cos(\omega_0 t)$ through the coefficient samples of the function and its first two derivatives. For OSNR Tier 3, the crucial change is semantic: Dirichlet, Neumann, and curvature boundary data become direct coefficient assignments. A clamped boundary at knot $k_b$, for example, is imposed by setting \[ \mathbf{c}[k_b] = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \] rather than by adding a soft loss term. The same vector-valued structure also produces the block-circulant Hermite Gram system derived in Appendix~\ref{app:hermite-block-gram}. \section{OSNR Architecture} \subsection{Master core apparatus} The current OSNR research code separates into three production tiers. Each tier has a different admissible function space, solver topology, and verification target. \begin{center} \scriptsize \begin{tabular}{@{}p{0.3\linewidth}p{0.3\linewidth}p{0.3\linewidth}@{}} \toprule \multicolumn{3}{c}{\textbf{OSNR Master Packaging Core}} \\ \midrule \textbf{Tier 1: Steady-State} & \textbf{Tier 2: Adaptive} & \textbf{Tier 3: Neural Operator} \\ \midrule Uniform E-splines & Non-uniform shift-splines & 2D tensor-product Hermite \\ $O(M\log M)$ circulant FFT & TLS matrix-pencil FRI & 9-stream local calculus \\ Fixed pole-locking inversion & Oblique cross-Gram shield & Parallel $9\times9$ DFT solver \\ \midrule Profile: $314.86$ dB, $4.21$ ms & Profile: $118.78$ dB, $97.7\%$ sparsity & Profile: $183.77$ dB, $1.28$ ms 2D CFD \\ \bottomrule \end{tabular} \end{center} \subsection{Tier 1: deterministic operator-bound systems} Tier 1 targets smooth physical fields that lie in, or close to, the null space of a known linear operator. Examples include Helmholtz waves, damped oscillators, and linear constant-coefficient PDE components. The solver pipeline is: \begin{enumerate}[leftmargin=2em] \item identify $L$ and its poles $\bs{\alpha}$; \item construct calibrated E-spline or null-space dictionaries; \item recover coefficients through a stable linear solve or circulant FFT inversion; \item compute differential quantities through operator identities rather than autograd. \end{enumerate} \subsection{Tier 2: adaptive sparse continua} Tier 2 targets composite fields with sparse discontinuities or shocks: \[ s = s_{\mathrm{smooth}} + s_{\mathrm{sparse}}. \] Uniform grids are not sufficient for non-bandlimited discontinuities at sub-grid coordinates. The sparse tier must first identify finite-rate innovations, then adapt the dictionary to those coordinates: \begin{enumerate}[leftmargin=2em] \item estimate innovation locations using FRI or matrix-pencil methods; \item snap sparse knots to the recovered coordinates; \item solve a cross-Gram-coupled sparse-plus-smooth inverse problem; \item debias with scale-invariant ridge stabilization. \end{enumerate} \subsection{Tier 3: higher-order neural operators} Tier 3 targets families of PDE solutions rather than a single fitted field. Neural operators such as DeepONet and MIONet learn maps between function spaces by pairing branch networks, which encode input functions or boundary data, with trunk networks, which encode query coordinates \cite{lu2021deeponet,jin2022mionet}. Spline-PINN shows a related but distinct path: a CNN predicts Hermite spline coefficients, and a continuous Hermite spline layer evaluates PDE residuals without finite-difference losses \cite{wandel2022splinepinn}. OSNR adopts the continuous Hermite idea but removes the black-box coordinate trunk where the governing operator and boundary calculus are known. Let \[ \Phi(t) = \begin{bmatrix} \phi_0(t) & \phi_1(t) & \phi_2(t) \end{bmatrix}^{\top} \] be a second-order Hermite generator compactly supported on $[-1,1]$. The generalized higher-order Hermite construction of Schmitter, Badoual, Uhlmann, Fageot, and Unser stores value, slope, and curvature data at each knot and interpolates all three channels exactly \cite{schmitter2016hermite}: \[ \phi_0^{(r)}(k)=\delta_{r0}\delta_k,\qquad \phi_1^{(r)}(k)=\delta_{r1}\delta_k,\qquad \phi_2^{(r)}(k)=\delta_{r2}\delta_k, \qquad r=0,1,2. \] The synthesized field is \[ f(t) = \sum_{k\in\Z} c_0[k]\phi_0(t-k) + c_1[k]\phi_1(t-k) + c_2[k]\phi_2(t-k), \] where \[ \mathbf{c}[k] = \begin{bmatrix} c_0[k] \\ c_1[k] \\ c_2[k] \end{bmatrix} = \begin{bmatrix} f(k) \\ f'(k) \\ f''(k) \end{bmatrix} \] for exactly interpolated data. The Hermite generator in \cite{schmitter2016hermite} is piecewise polynomial-exponential, $C^2$, compactly supported, and reproduces cubic polynomials as well as trigonometric functions. This is the missing boundary mechanism for a neural-operator OSNR tier: Dirichlet, Neumann, and curvature constraints are coefficient assignments, not soft loss penalties. For a PDE solution operator \[ \mathcal{S}: u \mapsto v, \] the branch-side network or analytic encoder should output Hermite coefficient tensors \[ \{\mathbf{c}[k;u]\}_{k\in\Z}, \] while the trunk side is replaced by deterministic Hermite synthesis. In multidimensional domains, tensor products of one-dimensional Hermite generators give mixed value/derivative channels, matching the construction used by Spline-PINN for continuous PDE residuals \cite{wandel2022splinepinn}. Unlike Spline-PINN, OSNR evaluates continuous energies and cross-correlations in coefficient space through the block-Gram calculus derived in Appendix~\ref{app:hermite-block-gram}, avoiding Monte Carlo residual quadrature whenever the operator and boundary model admit exact inner products. \section{Continuous-Discrete Calibration} The most important implementation condition is the cardinal coordinate map, which preserves the shift-invariant structure required by spline filtering calculus \cite{unser2005cardinal1,unser2005cardinal2}. \[ v_k(x) = \frac{x}{T} - k. \] Here $T$ is the physical knot spacing. For a domain $[0,D]$ and a compact generator of support order $q$, a stable finite dictionary uses \[ T = \frac{D}{M-q}. \] Then \[ v_k(x+T) = v_k(x)+1, \] so one physical knot step maps to one cardinal interval. \begin{proposition}[Failure of uncalibrated grid refinement] If the map is implemented as $v_k(x)=x+b_k$ with biases distributed over a fixed interval while $M$ changes, increasing $M$ decreases the relative shift between adjacent columns without changing physical support. The resulting dictionary columns become nearly collinear, and the Gram matrix $\A^\top \A$ develops near-zero eigenvalues. \end{proposition} \begin{proof}[Sketch] Let adjacent atoms be sampled as $\phi(x+b_k)$ and $\phi(x+b_{k+1})$. If $b_{k+1}-b_k=O(1/M)$ while the support width of $\phi$ is fixed, then a first-order expansion gives \[ \phi(x+b_{k+1}) = \phi(x+b_k) + O(1/M). \] Thus adjacent columns converge to each other as $M$ increases. The Gram matrix approaches rank deficiency and the pseudoinverse amplifies roundoff along small singular directions. \end{proof} This explains why a derivative residual can be zero while reconstruction fails. If the derivative dictionary is defined algebraically by $\A_{d2}=-k^2 \A$, then the PDE residual \[ \A_{d2}\cvec + k^2 \A \cvec \] is identically zero regardless of whether $\A$ is a well-conditioned reconstruction basis. \subsection{2D tensor-product Hermite expansion} The 2D Tier 3 engine is obtained by tensorizing the one-dimensional second-order Hermite streams. Let \[ h_i^x(x),\qquad h_j^y(y),\qquad i,j\in\{0,1,2\}, \] denote the value, slope, and curvature Hermite generators along the two axes. The tensor-product basis functions are \[ h_{i,j}(x,y) = h_i^x(x)h_j^y(y), \qquad i,j\in\{0,1,2\}. \] For a grid node $(k,\ell)$, OSNR stores a nine-stream local state \[ \mathbf{c}_{k,\ell} = \begin{bmatrix} f & \partial_x f & \partial_{xx}f & \partial_y f & \partial_{xy}f & \partial_{xxy}f & \partial_{yy}f & \partial_{xyy}f & \partial_{xxyy}f \end{bmatrix}_{(k,\ell)}^{\top}. \] The synthesized field is \[ f(x,y) = \sum_{k,\ell} \sum_{i=0}^{2}\sum_{j=0}^{2} c_{i,j}[k,\ell]\, h_i^x\!\left(\frac{x}{T_x}-k\right) h_j^y\!\left(\frac{y}{T_y}-\ell\right). \] This formula is the tensor-product analogue of the Schmitter et al. Hermite generator and the multidimensional counterpart of the periodic inner-product calculus of Badoual, Schmitter, and Unser. The practical consequence is that the spatial calculus of 2D physical fields is a forward coefficient operation. In the stream-function formulation for incompressible flow, \[ v_x=\partial_y a_z,\qquad v_y=-\partial_x a_z, \] and therefore \[ \nabla\cdot \mathbf{v} = \partial_x\partial_y a_z-\partial_y\partial_x a_z = 0 \] up to the commutation error of the discrete finite-difference stencil. The nonlinear transport term is evaluated as \[ (\mathbf{v}\cdot\nabla)\mathbf{v} = \begin{bmatrix} v_x\partial_x v_x+v_y\partial_y v_x\\ v_x\partial_x v_y+v_y\partial_y v_y \end{bmatrix}, \] and viscous diffusion as \[ \Delta\mathbf{v} = \begin{bmatrix} \partial_{xx}v_x+\partial_{yy}v_x\\ \partial_{xx}v_y+\partial_{yy}v_y \end{bmatrix}. \] All operators appearing in these expressions are evaluated by shift-invariant finite-difference ladders on the Hermite coefficient streams during the forward pass. No backward-mode automatic differentiation tape is constructed; the measured PyTorch autograd graph allocation in all Tier 3 validations is therefore $0.00$ bytes. \section{Autograd-Free Differential Calculus} \label{sec:autograd-free} For an exponential spline with pole vector $\bs{\alpha}$, applying a first-order operator $(D-\alpha_m)$ reduces the order of the spline \cite{unser2005cardinal1,delgadogonzalo2012exponential}: \[ (D-\alpha_m)\beta_{\bs{\alpha}}(t) = \beta_{\bs{\alpha}\setminus \alpha_m}(t) - e^{\alpha_m} \beta_{\bs{\alpha}\setminus \alpha_m}(t-1). \] This is the finite-difference ladder. Differential fields can be evaluated by filtering coefficients or by applying deterministic lower-order dictionary maps. For the Helmholtz null space, \[ \frac{\dd^2}{\dd x^2} s(x) = -k^2 s(x) \] inside the smooth spans, with boundary innovations handled separately. This removes the need for backward-mode automatic differentiation in Tier 1 PDE residual evaluation. The differential operator is encoded in the basis and coefficient algebra. \section{Circulant and FFT Solvers} When the dictionary is shift-invariant and periodized, the Gram matrix is circulant. This is the coefficient-domain form of the exact periodic inner-product calculus for spline curves and functions \cite{badoual2016inner,badoual2018periodic}: \[ \G = \begin{bmatrix} g_0 & g_{M-1} & \cdots & g_1 \\ g_1 & g_0 & \cdots & g_2 \\ \vdots & \vdots & \ddots & \vdots \\ g_{M-1} & g_{M-2} & \cdots & g_0 \end{bmatrix}. \] Such matrices are diagonalized by the discrete Fourier transform: \[ \G = \mathbf{F}^{-1} \Lambda \mathbf{F}. \] Solving $\G \cvec = \mathbf{b}$ reduces to \[ \widehat{\cvec}[\ell] = \frac{\widehat{\mathbf{b}}[\ell]}{\lambda_\ell}. \] This converts $O(M^3)$ dense solves into $O(M\log M)$ FFT operations. \section{Tomographic Radiance Fields as Spline Inverse Problems} The failed direct ray-kernel DL3DV experiment clarifies an important modeling boundary. A nontrivial view-synthesis scene is not naturally a smooth map from ray origin and direction to RGB. The physically shared quantity is a latent field in space, observed through line or ray measurements. The spline tomography literature gives the appropriate replacement model: represent the unknown continuous field by shifted basis functions, push those basis functions through the forward projector, and solve the coefficient inverse problem with fast adjoint and normal operators \cite{nilchian2013fast,mccann2016fast,donati2018multiscale,haouchat2025generalized}. Jin et al. make the complementary point that inverse problems whose normal operators are convolutional admit physics-aware direct inversions before any learned artifact-removal stage \cite{jin2017deep}; in OSNR, the direct inverse is the primary object, and any neural residual must remain secondary. For a first linearized density or opacity stage, write \[ \sigma(\mathbf{x})=\sum_{\mathbf{k}} a[\mathbf{k}]\,\varphi_\sigma(\mathbf{x}-\Lambda\mathbf{k}), \qquad \mathbf{g}=H\mathbf{a}+\boldsymbol{\eta}, \] where $H$ samples line or ray integrals of the spline density field. The regularized inverse update is \[ \mathbf{a}^{\star} = \arg\min_{\mathbf{a}} \frac12\|W^{1/2}(H\mathbf{a}-\mathbf{g})\|_2^2 +\lambda R(\mathbf{a}), \] with $W$ encoding ray confidence or frequency reliability, as in weighted phase-retrieval formulations \cite{bostan2016variational}. For quadratic $R(\mathbf{a})=\|L\mathbf{a}\|_2^2$, the normal equation is \[ (H^\top W H+\lambda L^\top L)\mathbf{a} = H^\top W\mathbf{g}. \] The computational opportunity is the McCann--Donati normal-operator identity. For shift-invariant basis functions and locally stationary projection blocks, $H^\top H$ acts as a discrete convolution: \[ (H^\top H\mathbf{a})[\mathbf{k}] = (\mathbf{a}*\mathbf{r})[\mathbf{k}], \] where $\mathbf{r}$ is the sampled autocorrelation of the projected basis function. Thus the expensive repeated normal-operator application inside conjugate gradients or ADMM becomes a Fourier-domain multiplication. Multiscale basis functions then provide a controlled coarse-to-fine path that is robust to pose and angular uncertainty \cite{donati2018multiscale}. Fast forward projection and sparse acquisition variants can be imported from the same lineage: Arcadu et al. use Fourier regridding with minimal oversampling for efficient forward projectors, while Donati et al. show how randomized STEM sampling can be coupled to regularized tomographic recovery \cite{arcadu2016forward,donati2017compressed}. This does not make full NeRF rendering linear. The volume-rendering equation contains transmittance \[ C(r)=\int T(t)\sigma(r(t))c(r(t),\mathbf{d})\,dt, \qquad T(t)=\exp\!\left(-\int_0^t\sigma(r(s))\,ds\right), \] so exact RGB fitting remains nonlinear in $\sigma$. The OSNR route is therefore staged: first recover coarse density/support with a tomographic spline inverse solve; then refine the density multiscale; then solve color/radiance coefficients on the recovered support; and finally apply visibility-weighted nonlinear corrections. Positivity, support constraints, total-variation, Hessian, and sparse-innovation priors enter naturally through the constrained ADMM machinery developed for spline tomography \cite{nilchian2013constrained,nilchian2015spline}. The controlled runner \texttt{apps\_industrial\_breakthrough/spline\_tomographic\_radiance\_solver.py} validates only the linearized operator claim. It builds a periodic synthetic density field, samples $48$ discrete projection directions, constructs $H^\top H$ explicitly once from a delta impulse, and then replaces all subsequent normal-operator applications by FFT convolution. On a $96\times96$ field, the convolutional normal operator matches explicit $H^\top H$ with relative error $2.4462\times10^{-7}$. One explicit normal-operator application costs $2.4956$ ms, while the FFT version costs $0.0969$ ms. Solving the same ridge-regularized inverse problem by conjugate gradients takes $188.1746$ ms with explicit normals and $4.9666$ ms with FFT normals, yielding a reconstruction PSNR of $22.3578$ dB under noisy sparse projections. This is not yet a NeRF result; it is the isolated mathematical validation that the spline-tomographic normal operator can be diagonalized as the literature predicts. The next controlled runner, \texttt{apps\_industrial\_breakthrough/spline\_ray\_operator\_validation.py}, validates the more fundamental Haouchat-style requirement: the forward ray operator and its adjoint must be matched before any real-scene radiance experiment is meaningful. On a $36\times36$ coefficient grid with $2688$ parallel rays, the script constructs two explicit small operators for auditability: a pixel basis and a quadratic tensor-product spline basis. The spline data are generated by the spline operator itself, and both models solve the same noisy inverse problem with conjugate gradients. The adjoint identity $\langle Hc,p\rangle=\langle c,H^\top p\rangle$ holds to $1.5672\times10^{-15}$ relative error for the spline operator and $2.7427\times10^{-16}$ for the pixel operator. At the same coefficient count, the spline inverse reconstructs the continuous rendered target at $54.0869$ dB, while the pixel model reaches only $27.0125$ dB. This is an intentionally controlled operator test: it proves that the coefficient-domain ray basis and adjoint are now correctly formulated, not that the full DL3DV visibility problem is solved. The basis choice itself is not incidental. Following the exponential-spline construction of Delgado-Gonzalo, Thevenaz, and Unser \cite{delgadogonzalo2012exponential}, a one-dimensional cardinal exponential B-spline associated with poles $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_N)$ has Fourier-domain form \[ \widehat{\beta}_{\boldsymbol{\alpha}}(\omega) = \prod_{m=1}^{N} \frac{1-\exp(\alpha_m-j\omega)}{j\omega-\alpha_m}. \] The two-dimensional smooth tier then uses the tensor-product generator \[ \varphi_{\boldsymbol{\alpha}_x,\boldsymbol{\alpha}_y}(x,y) = \beta_{\boldsymbol{\alpha}_x}(x)\, \beta_{\boldsymbol{\alpha}_y}(y), \] so the basis can reproduce the local modes implied by the operator rather than merely interpolate samples. The controlled sweep \texttt{apps\_industrial\_breakthrough/exponential\_spline\_basis\_sweep.py} tests this precision-first hypothesis on the same $1920$ rays and $784$ coefficients while changing only the tensor-product basis. The data are generated by a damped-harmonic exponential spline with poles $(-\lambda,-\lambda+j\omega,-\lambda-j\omega)$, and all candidate operators reuse the identical ray geometry. The matched damped-harmonic basis reaches $56.3837$ dB, compared with $52.1420$ dB for a monotone exponential-decay basis, $50.0620$ dB for an undamped harmonic basis, $49.9084$ dB for a quadratic polynomial spline, and $28.7579$ dB for a pixel box basis. This isolates the important design rule for the next radiance-field stage: the highest precision should come from tensor-product operator splines whose poles are matched to the expected local dynamics, with sparsity and compression applied after the operator basis is correct. The follow-up runner \texttt{apps\_industrial\_breakthrough/exponential\_spline\_pole\_sweep.py} turns this from a hand-picked basis comparison into a deterministic pole-selection problem. It fixes the target field, ray geometry, coefficient count, noise level, ridge parameter, and conjugate-gradient budget, then sweeps damped-harmonic exponential splines of orders $2$, $3$, and $4$ over $\lambda\in\{0.18,0.28,0.35,0.42,0.56,0.72\}$ and periods $\{8,10,12,16\}$. The target is generated by the order-$3$ pole set $(-0.42,-0.42+j2\pi/10,-0.42-j2\pi/10)$. The pole sweep correctly ranks that matched operator first at $55.2077$ dB on $1280$ rays and $576$ coefficients. The best order-$4$ candidate reaches $52.4859$ dB, and the best order-$2$ candidate reaches only $42.1900$ dB. This is the first automated evidence that pole selection is a meaningful OSNR model-selection axis: increasing support/order blindly does not dominate, while matching the operator poles controls reconstruction precision. The support/regularity runner \texttt{apps\_industrial\_breakthrough/exponential\_spline\_support\_regularization\_sweep.py} then isolates the opposite regime: a field generated by the shortest first-order Green-matched spline with pole $\alpha=-0.42$. The one-pole basis has support length $1$ and no continuity guarantee, but it exactly matches the local Green mode. It reaches $65.7918$ dB with operator density $0.0362$ and an average of $20.86$ active coefficients per ray. The smooth order-$3$ mixed basis $[0,0,\alpha]$ reaches only $27.1299$ dB and has density $0.1084$ with $62.44$ active coefficients per ray. This confirms a second design rule: if the modeled object is a Green response or sparse innovation, the shortest matched spline can be both more accurate and more localized than a smoother high-order basis. Regularity should be introduced because the signal class requires it, not by default. Finally, \texttt{apps\_industrial\_breakthrough/exponential\_spline\_basis\_selection\_map.py} evaluates the full pole-multiset selection problem. It tests seven target regimes against fourteen candidate bases: first-order Green response, repeated real poles, two distinct real poles, zero-augmented smooth operators, damped oscillators, pure polynomial splines, and support-$4$ mixed bases. Across all regimes, the exact pole multiset ranks first. This is the strongest evidence so far that OSNR basis design should be formulated as operator pole selection rather than degree selection. Order controls support and regularity; pole multiplicity and location control the reproduced null-space modes. Higher order improves asymptotic approximation power for smooth functions, but it is not a substitute for matching the operator that generated the signal. The final synthetic step, \texttt{apps\_industrial\_breakthrough/exponential\_spline\_operator\_inference.py}, removes access to rendered target PSNR during model selection. Each candidate basis is fitted on $75\%$ of the rays and scored on held-out rays using a normalized validation residual plus small support and density penalties. This exposes a practical distinction between generative pole matching and predictive operator selection. Under sparse-ray training, the oracle PSNR basis is sometimes a smoother support-$4$ model rather than the exact generating basis, because the added regularity improves interpolation across unobserved rays. The validation score still selects the oracle basis in four of seven regimes and stays within $0.2039$ dB of oracle in all cases, with mean PSNR loss $0.0465$ dB. Thus the operational rule becomes: use the differential operator poles as the first prior, then choose among nearby pole augmentations by held-out measurement prediction rather than training residual. We then stress-test the same idea in a mixed local-operator field with \texttt{apps\_industrial\_breakthrough/exponential\_spline\_local\_operator\_adaptation.py}. The synthetic field assigns different pole multisets to different spatial regions and compares a single global basis, a held-out-ray local block selector, and an oracle local label map. The oracle local solve reaches $33.5291$ dB, while the best global held-out model reaches $28.6559$ dB. This proves that local operator adaptation has substantial headroom. However, a blind ray-only greedy block selector reaches only $28.8074$ dB with $25\%$ block-label accuracy. Global line measurements make small independent block labels weakly identifiable unless the selection objective includes stronger spatial priors, localized measurements, or joint segmentation/coefficient optimization. The sparse-innovation remedy is tested in \texttt{apps\_industrial\_breakthrough/exponential\_spline\_fri\_region\_adaptation.py}. Instead of selecting arbitrary blocks, the method first reconstructs a global proxy, computes derivative-energy profiles, extracts sparse FRI-style transition proposals, expands them into a small boundary lattice, and then jointly scores region geometry and region pole choices on held-out rays. The raw derivative peaks locate approximate boundaries at $x=(-1.4149,1.4149)$ and $y=3.0319$; held-out refinement moves them to $x=(-3.4149,3.4149)$ and $y=4.0319$, yielding $100\%$ region-label agreement at the coefficient-grid resolution. The resulting FRI-region adaptive model reaches $34.2000$ dB, outperforming both the global model and the nominal oracle-region pole assignment. This is the first successful local adaptation mechanism: sparse innovation proposals supply the missing spatial prior that held-out global rays alone could not provide. The robustness sweep \texttt{apps\_industrial\_breakthrough/exponential\_spline\_fri\_region\_stress\_sweep.py} repeats the experiment over three projection-angle budgets and three additive noise levels. To keep the sweep diagnostic rather than combinatorial, it uses a compact region-basis candidate set containing the physical oracle family, the single-case selected family, a smoother support-$4$ family, and the best global family; the exhaustive $5^4$ region-basis search remains available as an optional mode. Across all nine stress cases, the FRI-region model improves over the best global held-out basis. The gain increases with measurement density, from a mean $+1.1667$ dB at $12$ angles to $+5.0277$ dB at $24$ angles, while the recovered region accuracy rises from $89.50\%$ to $100.00\%$. This confirms that the sparse-innovation step is not a one-off artifact: as the inverse problem receives enough projections to identify the transition set, local operator adaptation becomes reliably beneficial. We then deliberately break the axis-aligned assumption with \texttt{apps\_industrial\_breakthrough/exponential\_spline\_fri\_curve\_adaptation.py}. The target field is generated by three smooth transition curves: two vertical pole-boundary curves and one top-interface curve. A row/column FRI-style derivative tracker fits low-order curve proposals from the global proxy and refines them by held-out rays. This harder test exposes the current bottleneck. The global held-out basis reaches $32.2114$ dB, while the true curved local operator assignment reaches $41.5032$ dB, proving that curved local operators have large headroom. However, the detected curved selector reaches only $31.1969$ dB despite $91.00\%$ region-label agreement and boundary RMSE $1.4606$. The failure is not the absence of local operator advantage; it is the scoring layer. Under curved imperfect labels, the held-out ray residual prefers smoother surrogate pole assignments rather than the physical local pole map. The next algorithmic step is therefore joint geometry--basis--coefficient refinement, or a region-contrastive validation score that prevents the local operator assignment from collapsing to a globally smooth surrogate. The joint refinement runner \texttt{apps\_industrial\_breakthrough/exponential\_spline\_fri\_curve\_joint\_refinement.py} implements this next correction. It expands the curve-offset lattice near the best FRI proposal, fits coefficients for each candidate local operator assignment, and augments the held-out residual with an edge-consistency contrast term that rewards reconstructions whose gradient energy concentrates on the proposed sparse transition curves. This converts the curved selector from a negative result into a partial recovery: the selected joint model reaches $34.7216$ dB, a $+2.5102$ dB gain over the global held-out basis. The best candidate present in the searched family reaches $35.8310$ dB, while the true-curve oracle remains at $41.5032$ dB. Thus the scoring fix is directionally correct but incomplete. The remaining gap now separates two effects: curve localization error and the limited pole-assignment candidate family. This is the cleanest current target for further algorithmic work. Finally, \texttt{apps\_industrial\_breakthrough/exponential\_spline\_geometry\_reproduction.py} isolates the more geometric point raised by the exponential-spline curve and surface literature \cite{delgadogonzalo2012exponential}. The target is a closed harmonic curve with modes up to order three, represented from only twelve parameter samples and eight control points. The matched compact harmonic E-spline uses the pole set $\{0,\pm j2\pi/M,\pm j4\pi/M,\pm j6\pi/M\}$ and reaches dense curve RMSE $5.4457\times10^{-6}$. With the same number of control points and the same samples, a generic cubic polynomial spline reaches only $1.0256\times10^{-2}$ RMSE, and a piecewise-linear polygon reaches $3.8693\times10^{-2}$ RMSE. This is a small but important result: if the expected geometry is known to be harmonic, elliptic, spherical, cylindrical, or otherwise parametrizable by a known exponential-polynomial family, then OSNR should place that family directly in the geometric span rather than recover it indirectly through a generic volumetric grid. For NeRF-like scenes this suggests a patch-based route: segment or initialize object surfaces with geometry-reproducing parametric E-splines, then fit texture/radiance on those surfaces. For PDE domains it suggests an even cleaner route: represent both boundary geometry and the solution field in operator-matched spline spaces. The follow-up optimizer \texttt{apps\_industrial\_breakthrough/eggroll\_spline\_shape\_optimizer.py} tests whether this geometric advantage can be used when the target shape is not known in closed form. Motivated by Schmitter and Unser's continuous-domain shape projectors and functional PCA construction \cite{schmitter2018landmark}, the experiment represents a closed spline curve by $16$ control points but restricts the learned geometric search to an $8$-dimensional continuous shape subspace. This is the analogue of replacing an arbitrary coordinate-field parameter vector by a learned spline-shape chart. The stochastic search component is motivated by the EGGROLL low-rank evolution-strategy result \cite{sarkar2026eggroll}: rank-one perturbations can be evaluated as hardware-friendly low-rank updates, but the experiment separates this hardware trick from the geometric prior itself. The result is deliberately diagnostic. Full Gaussian ES over all $32$ control coordinates reaches dense RMSE $2.9097\times10^{-2}$ after $67{,}200$ forward evaluations, while rank-one EGGROLL-style perturbations applied directly to the raw control matrix reach $3.3517\times10^{-2}$. Thus low-rank noise alone does not solve geometry discovery. When the same evaluation budget is spent inside the Schmitter-style spline subspace, the error drops to $6.9268\times10^{-3}$; rank-one EGGROLL perturbations inside that subspace reach a comparable $7.4862\times10^{-3}$. The algebraic continuous-subspace projection oracle reaches $3.5689\times10^{-4}$ in $0.293$ ms, exposing the remaining optimization gap. The practical implication is precise: the promising route is not blind evolution over arbitrary OSNR coefficients, but variable projection. Use stochastic low-rank search only for nonlinear geometry, visibility, and knot variables; solve the linear radiance or texture coefficients algebraically once a candidate geometry is proposed. The next controlled runner, \texttt{apps\_industrial\_breakthrough/eggroll\_adjoint\_variable\_projection.py}, implements that variable-projection step explicitly. A one-dimensional spline boundary partitions a $40\times40$ radiance field into two continuous regions. For each candidate boundary, the code constructs the ray operator $H(\varphi)$ by projecting masked smooth atoms, eliminates the linear radiance coefficients by the adjoint normal equation \[ c^\star(\varphi)=\left(H(\varphi)^\top H(\varphi)+\lambda I\right)^{-1}H(\varphi)^\top y, \] and scores the resulting field on acquisition directions not used in the coefficient solve. This is the minimal inverse-problem analogue of a NeRF geometry/radiance separation: nonlinear geometry is searched, while linear radiance is solved in closed form. The experiment clarifies both the opportunity and the bottleneck. The mean-geometry variable-projection baseline reaches $16.8968$ dB on held-out projections. Full Gaussian ES over raw boundary controls improves to $18.9082$ dB, rank-one EGGROLL over raw controls reaches $19.9034$ dB, and rank-one EGGROLL in the six-dimensional spline subspace reaches $20.1043$ dB. However, the true-geometry variable-projection oracle reaches $33.1874$ dB with a field RMSE of $3.3521\times10^{-2}$. Thus the adjoint variable-projection mechanism is working, but stochastic boundary discovery remains underidentified from the current projection residual alone. This is an important negative constraint for the NeRF/SIREN campaign: the next improvement must add stronger geometry evidence, such as FRI edge measurements, silhouette consistency, epipolar visibility constraints, or a learned continuous shape prior. More low-rank perturbation budget alone is unlikely to close the oracle gap. The Haouchat-matched follow-up \texttt{apps\_industrial\_breakthrough/haouchat\_matched\_variable\_projection.py} performs that correction. Instead of using a primitive projection mask, it builds the inner ray operator from the same quadrature-evaluated tensor-product spline basis used in the $52$--$56$ dB matched-ray experiments. The target data are generated by a damped-harmonic exponential spline ray operator. Candidate geometries still define a boundary-dependent masked coefficient dictionary, but the forward map is now \[ H_\varphi = H_{\beta_L}\,\Phi(\varphi), \] where $H_{\beta_L}$ is the Haouchat-style ray projector for the selected tensor-product basis and $\Phi(\varphi)$ applies the boundary-dependent coefficient atoms. The stochastic outer loop is also given a weak FRI-like edge observation of the boundary, so the score combines held-in projection residual and edge consistency. This is the first experiment in this line that combines all three ingredients: matched spline rays, adjoint variable projection, and sparse boundary evidence. The result changes the interpretation sharply. With the matched damped-harmonic operator, the true-geometry projection oracle reaches $96.9888$ dB on held-out rays, confirming that the inner ray/inverse model itself is not the limiting factor. The edge-aware rank-one spline-subspace EGGROLL search reaches $42.2339$ dB, up from the mean-geometry baseline of $33.3709$ dB, with boundary RMSE reduced from $9.5169\times10^{-2}$ to $2.7066\times10^{-2}$. The exponential-decay candidate also benefits from edge evidence, improving from $40.6125$ dB without the edge term to $41.0644$ dB with it. This supports the current thesis: OSNR does not need a dense NeRF-style MLP to represent the radiance once the operator is matched; the hard remaining problem is physically constrained geometry and visibility discovery. The follow-up \texttt{apps\_industrial\_breakthrough/haouchat\_fri\_edge\_variable\_projection.py} removes the remaining artificial part of the edge-aware score. Instead of injecting a boundary hint directly from the true geometry, it forms a measurement-derived FRI proxy: training rays are backprojected through the matched adjoint, row-wise derivatives of the normalized adjoint image are localized, and the resulting peak track is smoothed into a candidate boundary. This is not a complete multidimensional FRI surface solver, but it is a measurement-only sparse-transition proposal. In the default run, the synthetic edge hint has boundary RMSE $1.3676\times10^{-2}$, while the adjoint-derived edge track has RMSE $1.9117\times10^{-2}$. Using this measurement-derived edge evidence, the matched damped-harmonic model reaches $41.5632$ dB on held-out rays, compared with $33.3709$ dB for mean geometry and $96.9888$ dB for the true-geometry oracle. The cleaner synthetic edge hint reaches $45.9174$ dB in the same runner. The gap between $41.56$ dB and $45.92$ dB is useful: it quantifies the price of deriving geometry evidence from measurements rather than providing it externally. The experiment therefore validates the direction without hiding the remaining work. The next real-scene version should replace row-wise adjoint peaks by multi-view epipolar FRI proposals and visibility-aware surface clustering. The next controlled refinement, \texttt{apps\_industrial\_breakthrough/haouchat\_multiray\_fri\_edge\_projection.py}, replaces the row-wise adjoint peak estimate by a multi-ray residual selection loop. The adjoint boundary is used only as initialization. Each boundary control point is perturbed over a small local offset lattice, and candidates are scored by the held-in variable-projection residual plus curvature and anchor penalties: \[ \mathcal{S}(\varphi)= \|H_{\beta_L}\Phi(\varphi)c^\star(\varphi)-y_{\mathrm{score}}\|_2^2 +\eta\|\Delta^2\varphi\|_2^2 +\rho\|\varphi-\varphi_{\mathrm{adj}}\|_2^2. \] This converts the crude measurement edge into a ray-consistent FRI boundary proposal without accessing the hidden target geometry. The improvement is large. The row/adjoint edge has boundary RMSE $1.9117\times10^{-2}$; the multi-ray refined edge has RMSE $1.4848\times10^{-3}$, better than the noisy synthetic edge control ($1.3676\times10^{-2}$). With the matched damped-harmonic operator, direct variable projection on the multi-ray FRI boundary reaches $66.1488$ dB on held-out rays and field RMSE $2.5205\times10^{-3}$. The true-geometry oracle remains $96.9888$ dB, so the experiment is still controlled rather than a final SOTA benchmark, but it shows that the geometry bottleneck can be attacked algebraically by ray-consistent sparse innovation refinement. Interestingly, applying the stochastic EGGROLL search after this refined geometry is worse ($40.3169$ dB), so the current best path is not more random search but better deterministic boundary proposal. The robustness sweep \texttt{apps\_industrial\_breakthrough/haouchat\_multiray\_fri\_robustness\_sweep.py} repeats the deterministic part of the experiment across two boundary variants, projection-angle counts $\{12,18,24\}$, and additive noise levels $\{0,3\times10^{-4},10^{-3}\}$. The mean-geometry baseline averages $32.5922$ dB over the $18$ cases; direct row/adjoint edge projection averages $38.0604$ dB; multi-ray FRI variable projection averages $64.0770$ dB; and the true-geometry oracle averages $97.9883$ dB. The multi-ray method remains above $61.76$ dB in every tested case. This confirms that the $66.15$ dB result is not a one-off numerical accident, but a stable consequence of selecting boundary offsets by held-in ray residuals. \begin{table}[h] \centering \small \begin{tabular}{lcc} \toprule Quantity & Explicit normal & FFT convolution normal \\ \midrule Single $H^\top H$ application & $2.4956$ ms & $0.0969$ ms \\ CG inverse solve & $188.1746$ ms & $4.9666$ ms \\ Normal-operator relative error & \multicolumn{2}{c}{$2.4462\times10^{-7}$} \\ Reconstruction PSNR & \multicolumn{2}{c}{$22.3578$ dB} \\ \bottomrule \end{tabular} \caption{Controlled synthetic validation of the spline-tomographic normal-operator identity. The experiment isolates the linear inverse-problem component needed before returning to nonlinear radiance-field rendering.} \label{tab:spline-tomographic-radiance} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcc} \toprule Quantity & Pixel basis & Quadratic spline basis \\ \midrule Matched-adjoint relative error & $2.7427\times10^{-16}$ & $1.5672\times10^{-15}$ \\ CG inverse solve & $19.5137$ ms & $21.9203$ ms \\ Coefficient RMSE & $4.1936\times10^{-2}$ & $7.5915\times10^{-3}$ \\ Rendered reconstruction PSNR & $27.0125$ dB & $54.0869$ dB \\ \bottomrule \end{tabular} \caption{Controlled validation of the matched spline ray operator $H_\varphi$ and adjoint $H_\varphi^\top$ on $2688$ rays and $1296$ coefficients. The result establishes the correct operator foundation before porting the method to DL3DV camera rays with visibility weights.} \label{tab:spline-ray-operator-validation} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Basis profile & PSNR & Adjoint error & CG solve & Coefficient RMSE \\ \midrule Pixel box & $28.7579$ dB & $1.8997\times10^{-16}$ & $4.8182$ ms & $1.8724\times10^{-1}$ \\ Quadratic spline & $49.9084$ dB & $3.0959\times10^{-15}$ & $4.1478$ ms & $1.8928\times10^{-1}$ \\ Exponential decay & $52.1420$ dB & $2.0145\times10^{-16}$ & $4.7577$ ms & $1.1161\times10^{-2}$ \\ Harmonic exponential & $50.0620$ dB & $2.2573\times10^{-16}$ & $4.1359$ ms & $1.5208\times10^{-2}$ \\ Damped-harmonic exponential & $56.3837$ dB & $0.0000$ & $3.9541$ ms & $6.8398\times10^{-3}$ \\ \bottomrule \end{tabular} \caption{Tensor-product basis sweep for a matched ray inverse problem on $1920$ rays and $784$ coefficients. The target field is generated by the damped-harmonic exponential spline; all candidate bases reuse the same rays, regularization, and conjugate-gradient inverse solve.} \label{tab:exponential-spline-basis-sweep} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Best candidate class & Pole profile & PSNR & CG solve & Coefficient RMSE \\ \midrule Order $2$ & $\lambda=0.35$, period $16$ & $42.1900$ dB & $2.1103$ ms & $1.6017\times10^{-1}$ \\ Order $3$ & $\lambda=0.42$, period $10$ & $55.2077$ dB & $1.9365$ ms & $7.5642\times10^{-3}$ \\ Order $4$ & $\lambda=0.35$, period $16$ & $52.4859$ dB & $2.2178$ ms & $4.9702\times10^{-2}$ \\ \bottomrule \end{tabular} \caption{Deterministic pole-selection sweep for damped-harmonic tensor-product exponential splines. The target is generated by the order-$3$, $\lambda=0.42$, period-$10$ operator; the matched pole set ranks first across the tested grid.} \label{tab:exponential-spline-pole-sweep} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Pole multiset & Regularity & PSNR & CG solve & Density & Active/ray \\ \midrule $[\alpha]$ & $C^{-1}$ & $65.7918$ dB & $2.1925$ ms & $0.0362$ & $20.86$ \\ $[0]$ & $C^{-1}$ & $29.1917$ dB & $2.1154$ ms & $0.0362$ & $20.86$ \\ $[0,\alpha]$ & $C^0$ & $26.8630$ dB & $2.1370$ ms & $0.0724$ & $41.71$ \\ $[\alpha,\alpha]$ & $C^0$ & $26.5164$ dB & $2.1558$ ms & $0.0724$ & $41.71$ \\ $[0,0,\alpha]$ & $C^1$ & $27.1299$ dB & $2.1973$ ms & $0.1084$ & $62.44$ \\ $[\alpha,\alpha,\alpha]$ & $C^1$ & $26.8975$ dB & $2.0574$ ms & $0.1084$ & $62.44$ \\ \bottomrule \end{tabular} \caption{Support-versus-regularity sweep for a first-order Green target with $\alpha=-0.42$ on $1280$ rays and $576$ coefficients. The shortest matched basis wins because the target is Green-like rather than smooth.} \label{tab:exponential-spline-support-regularity} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Target regime & Best candidate & Best PSNR & Matched PSNR & Support & Regularity \\ \midrule Green $[\alpha]$ & $[\alpha]$ & $52.9697$ dB & $52.9697$ dB & $1$ & $C^{-1}$ \\ Repeated $[\alpha,\alpha]$ & $[\alpha,\alpha]$ & $45.6373$ dB & $45.6373$ dB & $2$ & $C^0$ \\ Two-real $[\alpha,\beta]$ & $[\alpha,\beta]$ & $44.6760$ dB & $44.6760$ dB & $2$ & $C^0$ \\ Smooth $[0,0,\alpha]$ & $[0,0,\alpha]$ & $49.4038$ dB & $49.4038$ dB & $3$ & $C^1$ \\ Damped oscillator & $[-\lambda,-\lambda\pm j\omega]$ & $48.4246$ dB & $48.4246$ dB & $3$ & $C^1$ \\ Polynomial $[0,0,0]$ & $[0,0,0]$ & $49.9511$ dB & $49.9511$ dB & $3$ & $C^1$ \\ Mixed $[0,0,\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $52.2576$ dB & $52.2576$ dB & $4$ & $C^2$ \\ \bottomrule \end{tabular} \caption{Pole-multiset basis-selection map on $768$ rays and $400$ coefficients. The exact pole multiset ranks first in every tested target regime, confirming that operator matching and pole multiplicity are distinct from simply increasing spline order.} \label{tab:exponential-spline-basis-selection-map} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Target regime & Oracle basis & Inferred basis & Oracle PSNR & Inferred PSNR \\ \midrule Green $[\alpha]$ & $[\alpha]$ & $[\alpha]$ & $35.3689$ dB & $35.3689$ dB \\ Repeated $[\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $[0,0,0,0]$ & $33.1575$ dB & $33.1202$ dB \\ Two-real $[\alpha,\beta]$ & $[0,0,\alpha,\alpha]$ & $[0,0,0,0]$ & $32.2104$ dB & $32.0064$ dB \\ Smooth $[0,0,\alpha]$ & $[0,0,\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $34.9124$ dB & $34.9124$ dB \\ Damped oscillator & $[0,0,\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $34.6249$ dB & $34.6249$ dB \\ Polynomial $[0,0,0]$ & $[0,0,0,0]$ & $[0,0,\alpha,\alpha]$ & $34.9254$ dB & $34.8408$ dB \\ Mixed $[0,0,\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $[0,0,\alpha,\alpha]$ & $36.1621$ dB & $36.1621$ dB \\ \bottomrule \end{tabular} \caption{Measurement-driven operator inference with a held-out ray split. The selected basis is chosen without rendered target access; it remains within $0.2039$ dB of the oracle PSNR basis across all tested regimes.} \label{tab:exponential-spline-operator-inference} \end{table} \begin{table}[h] \centering \small \begin{tabular}{@{}p{0.22\linewidth}p{0.34\linewidth}p{0.27\linewidth}p{0.10\linewidth}@{}} \toprule Model & PSNR & Selection signal & Block-label accuracy \\ \midrule Global held-out basis & $28.6559$ dB & Held-out rays & n/a \\ Local greedy basis & $28.8074$ dB & Held-out rays & $25.00\%$ \\ Oracle local labels & $33.5291$ dB & Ground-truth region labels & $100.00\%$ \\ \bottomrule \end{tabular} \caption{Local operator adaptation stress test on a mixed pole-multiset field. The oracle gap confirms that local bases can matter, while the weak greedy-label recovery identifies the next algorithmic bottleneck.} \label{tab:exponential-spline-local-operator-adaptation} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule Model & Boundary source & Region accuracy & PSNR \\ \midrule Global held-out basis & none & n/a & $28.6559$ dB \\ Blind greedy blocks & fixed grid & $25.00\%$ & $28.8074$ dB \\ Oracle region labels & ground truth & $100.00\%$ & $33.5291$ dB \\ FRI-region adaptation & derivative peaks + held-out refinement & $100.00\%$ & $34.2000$ dB \\ \bottomrule \end{tabular} \caption{FRI-guided local operator adaptation. Sparse-innovation boundary proposals convert the weak blockwise selection problem into a region-level operator-selection problem and recover the local-basis advantage.} \label{tab:exponential-spline-fri-region-adaptation} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Projection angles & Mean global PSNR & Mean FRI PSNR & Mean gain & Region accuracy \\ \midrule $12$ & $28.2410$ dB & $29.4078$ dB & $+1.1667$ dB & $89.50\%$ \\ $18$ & $28.7692$ dB & $31.9501$ dB & $+3.1809$ dB & $86.00\%$ \\ $24$ & $28.9903$ dB & $34.0180$ dB & $+5.0277$ dB & $100.00\%$ \\ \bottomrule \end{tabular} \caption{FRI-region stress sweep averaged over additive noise levels $\{0,5\times10^{-4},2\times10^{-3}\}$. The sparse-transition prior becomes more valuable as the ray geometry provides enough measurements to localize region boundaries.} \label{tab:exponential-spline-fri-region-stress-sweep} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule Model & Boundary model & Region accuracy & PSNR \\ \midrule Global held-out basis & none & n/a & $32.2114$ dB \\ FRI curved selector & fitted curves & $91.00\%$ & $31.1969$ dB \\ Oracle curved labels & ground-truth curves & $100.00\%$ & $41.5032$ dB \\ \bottomrule \end{tabular} \caption{Curved-interface local operator adaptation. The oracle gap confirms large local-operator headroom, while the selected model identifies the next bottleneck: held-out ray residuals alone are not sufficient to choose physical pole assignments under imperfect curved segmentation.} \label{tab:exponential-spline-fri-curve-adaptation} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule Model & Selection signal & Region accuracy & PSNR \\ \midrule Global held-out basis & held-out rays & n/a & $32.2114$ dB \\ FRI curved selector & held-out rays & $91.00\%$ & $31.1969$ dB \\ Joint curve refinement & residual + edge contrast & $92.25\%$ & $34.7216$ dB \\ Best searched candidate & hidden PSNR oracle & n/a & $35.8310$ dB \\ Oracle curved labels & ground-truth curves & $100.00\%$ & $41.5032$ dB \\ \bottomrule \end{tabular} \caption{Joint curved-interface refinement. Adding edge-consistency contrast to the measurement score recovers a useful local-operator gain, but the remaining oracle gap shows that curved sparse-innovation geometry and pole assignment still need joint refinement.} \label{tab:exponential-spline-fri-curve-joint-refinement} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule Geometry basis & Parameters & Dense RMSE & Max error \\ \midrule Matched harmonic E-spline, $L=3$ & $16$ & $5.4457\times10^{-6}$ & $9.2024\times10^{-6}$ \\ Generic cubic polynomial spline & $16$ & $1.0256\times10^{-2}$ & $1.9079\times10^{-2}$ \\ Piecewise-linear polygon & $16$ & $3.8693\times10^{-2}$ & $1.6177\times10^{-1}$ \\ \bottomrule \end{tabular} \caption{Geometry-reproduction benchmark for a closed harmonic curve from twelve parameter samples. The matched exponential-spline pole set places the target geometry in the span; generic polynomial and polygonal bases require more parameters to reach the same precision.} \label{tab:exponential-spline-geometry-reproduction} \end{table} \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Optimizer & Search dimension & Dense RMSE & Runtime & Structural compression \\ \midrule Full Gaussian ES on controls & $32$ & $2.9097\times10^{-2}$ & $46.708$ ms & $0.00\%$ \\ Rank-one EGGROLL on controls & $32$ & $3.3517\times10^{-2}$ & $39.477$ ms & $0.00\%$ \\ Schmitter spline-subspace ES & $8$ & $6.9268\times10^{-3}$ & $35.385$ ms & $75.00\%$ \\ Rank-one EGGROLL in spline subspace & $8$ & $7.4862\times10^{-3}$ & $35.693$ ms & $75.00\%$ \\ Continuous subspace projection oracle & $8$ & $3.5689\times10^{-4}$ & $0.293$ ms & $75.00\%$ \\ \bottomrule \end{tabular} \caption{Low-rank stochastic geometry search on a continuous spline shape family. The result separates the EGGROLL hardware mechanism from the Schmitter-style geometric prior: raw rank-one perturbations are not enough, while low-dimensional continuous spline shape coordinates produce the large error reduction.} \label{tab:eggroll-spline-shape-optimizer} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/eggroll_spline_shape_optimizer_outputs/eggroll_spline_shape_optimizer.png} \caption{Controlled spline-shape discovery benchmark. Red points are sparse noisy observations, gray curves mark the target where shown, and black curves show the recovered continuous spline shape for each optimizer.} \label{fig:eggroll-spline-shape-optimizer} \end{figure} \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Method & Search dim. & Held-out PSNR & Field RMSE & Boundary RMSE & Search time \\ \midrule Mean geometry + variable projection & $0$ & $16.8968$ dB & $2.2344\times10^{-1}$ & $1.4287\times10^{-1}$ & $0.00$ ms \\ Full Gaussian ES controls & $18$ & $18.9082$ dB & $1.7413\times10^{-1}$ & $3.9663\times10^{-1}$ & $7132.59$ ms \\ Rank-one EGGROLL controls & $18$ & $19.9034$ dB & $1.5539\times10^{-1}$ & $3.4238\times10^{-1}$ & $7333.54$ ms \\ Rank-one EGGROLL spline subspace & $6$ & $20.1043$ dB & $1.7349\times10^{-1}$ & $3.2396\times10^{-1}$ & $7177.82$ ms \\ True geometry projection oracle & $6$ & $33.1874$ dB & $3.3521\times10^{-2}$ & $0.0000$ & $0.00$ ms \\ \bottomrule \end{tabular} \caption{Adjoint variable-projection geometry benchmark. Each candidate boundary defines $H(\varphi)$, radiance coefficients are eliminated by a ridge normal solve, and quality is measured on held-out projection directions. The oracle gap shows that coefficient elimination is not enough; physical boundary evidence must be strengthened.} \label{tab:eggroll-adjoint-variable-projection} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/eggroll_adjoint_variable_projection_outputs/eggroll_adjoint_variable_projection.png} \caption{Variable-projection radiance benchmark. Left panel is the target field; subsequent panels show recovered fields for mean geometry, full ES, rank-one control EGGROLL, rank-one spline-subspace EGGROLL, and true-geometry oracle. Red curves mark the recovered boundary.} \label{fig:eggroll-adjoint-variable-projection} \end{figure} \begin{table}[h] \centering \small \begin{tabular}{llccc} \toprule Basis & Method & Held-out PSNR & Field RMSE & Boundary RMSE \\ \midrule Polynomial quadratic & Mean geometry & $33.3249$ dB & $3.8306\times10^{-1}$ & $9.5169\times10^{-2}$ \\ Polynomial quadratic & Edge-aware EGGROLL & $41.7015$ dB & $3.8036\times10^{-1}$ & $2.3912\times10^{-2}$ \\ Polynomial quadratic & True-geometry oracle & $48.9766$ dB & $3.7926\times10^{-1}$ & $0.0000$ \\ Exponential decay & Mean geometry & $33.3627$ dB & $6.6277\times10^{-2}$ & $9.5169\times10^{-2}$ \\ Exponential decay & Edge-aware EGGROLL & $41.0644$ dB & $3.3052\times10^{-2}$ & $2.4452\times10^{-2}$ \\ Exponential decay & True-geometry oracle & $52.4398$ dB & $5.1297\times10^{-3}$ & $0.0000$ \\ Damped harmonic & Mean geometry & $33.3709$ dB & $6.6046\times10^{-2}$ & $9.5169\times10^{-2}$ \\ Damped harmonic & EGGROLL, no edge term & $39.9315$ dB & $4.0312\times10^{-2}$ & $5.0651\times10^{-2}$ \\ Damped harmonic & Edge-aware EGGROLL & $42.2339$ dB & $2.9005\times10^{-2}$ & $2.7066\times10^{-2}$ \\ Damped harmonic & True-geometry oracle & $96.9888$ dB & $3.2065\times10^{-5}$ & $0.0000$ \\ \bottomrule \end{tabular} \caption{Haouchat-matched variable projection. The inner operator is a quadrature-evaluated tensor-product spline ray projector, while the outer loop searches boundary geometry. Matching the damped-harmonic exponential basis restores the very high oracle ceiling, and adding sparse edge evidence moves the stochastic search into the $42$ dB held-out regime.} \label{tab:haouchat-matched-variable-projection} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/haouchat_matched_variable_projection_outputs/haouchat_matched_variable_projection.png} \caption{Haouchat-matched variable projection for the damped-harmonic basis. Panels show the target, mean geometry, no-edge EGGROLL, edge-aware EGGROLL, and true-geometry oracle.} \label{fig:haouchat-matched-variable-projection} \end{figure} \begin{table}[h] \centering \small \begin{tabular}{llccc} \toprule Basis & Method & Held-out PSNR & Field RMSE & Boundary RMSE \\ \midrule Polynomial quadratic & Measurement FRI edge & $41.2085$ dB & $3.8043\times10^{-1}$ & $2.5332\times10^{-2}$ \\ Polynomial quadratic & True-geometry oracle & $48.9766$ dB & $3.7926\times10^{-1}$ & $0.0000$ \\ Exponential decay & Measurement FRI edge & $40.5459$ dB & $3.4998\times10^{-2}$ & $4.8217\times10^{-2}$ \\ Exponential decay & True-geometry oracle & $52.4398$ dB & $5.1297\times10^{-3}$ & $0.0000$ \\ Damped harmonic & Mean geometry & $33.3709$ dB & $6.6046\times10^{-2}$ & $9.5169\times10^{-2}$ \\ Damped harmonic & Synthetic edge hint & $45.9174$ dB & $2.1096\times10^{-2}$ & $1.9775\times10^{-2}$ \\ Damped harmonic & Measurement FRI edge & $41.5632$ dB & $3.3397\times10^{-2}$ & $4.2715\times10^{-2}$ \\ Damped harmonic & True-geometry oracle & $96.9888$ dB & $3.2065\times10^{-5}$ & $0.0000$ \\ \bottomrule \end{tabular} \caption{Measurement-derived FRI edge variable projection. The sparse boundary proposal is estimated from adjoint backprojection and derivative peak localization, not from the hidden target geometry. It recovers most of the useful edge-aware gain but remains below the cleaner synthetic edge hint, identifying multi-view FRI geometry extraction as the next bottleneck.} \label{tab:haouchat-fri-edge-variable-projection} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/haouchat_fri_edge_variable_projection_outputs/haouchat_fri_edge_variable_projection.png} \caption{Measurement-derived FRI edge variable projection for the damped-harmonic basis. Panels show target, mean geometry, synthetic-edge EGGROLL, measurement-edge EGGROLL, and true-geometry oracle.} \label{fig:haouchat-fri-edge-variable-projection} \end{figure} \begin{table}[h] \centering \small \begin{tabular}{llccc} \toprule Basis & Method & Held-out PSNR & Field RMSE & Boundary RMSE \\ \midrule Polynomial quadratic & Multi-ray FRI variable projection & $48.6908$ dB & $3.7930\times10^{-1}$ & $1.4848\times10^{-3}$ \\ Polynomial quadratic & True-geometry oracle & $48.9766$ dB & $3.7926\times10^{-1}$ & $0.0000$ \\ Exponential decay & Multi-ray FRI variable projection & $51.9659$ dB & $5.6990\times10^{-3}$ & $1.4848\times10^{-3}$ \\ Exponential decay & True-geometry oracle & $52.4398$ dB & $5.1297\times10^{-3}$ & $0.0000$ \\ Damped harmonic & Mean geometry & $33.3709$ dB & $6.6046\times10^{-2}$ & $9.5169\times10^{-2}$ \\ Damped harmonic & Row/adjoint FRI edge & $37.7852$ dB & $4.7983\times10^{-2}$ & $6.8656\times10^{-2}$ \\ Damped harmonic & Multi-ray FRI variable projection & $66.1488$ dB & $2.5205\times10^{-3}$ & $1.4848\times10^{-3}$ \\ Damped harmonic & EGGROLL after multi-ray edge & $40.3169$ dB & $3.6051\times10^{-2}$ & $3.9707\times10^{-2}$ \\ Damped harmonic & True-geometry oracle & $96.9888$ dB & $3.2065\times10^{-5}$ & $0.0000$ \\ \bottomrule \end{tabular} \caption{Multi-ray FRI edge refinement. Local boundary offsets are selected by held-in ray residuals after variable projection. The direct refined geometry nearly closes the oracle gap for the polynomial and exponential-decay bases and raises the matched damped-harmonic profile to $66.15$ dB without synthetic edge injection.} \label{tab:haouchat-multiray-fri-edge-projection} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/haouchat_multiray_fri_edge_projection_outputs/haouchat_multiray_fri_edge_projection.png} \caption{Multi-ray FRI edge refinement for the damped-harmonic basis. The direct multi-ray FRI geometry, not the subsequent EGGROLL search, gives the dominant quality gain.} \label{fig:haouchat-multiray-fri-edge-projection} \end{figure} \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule Method & Mean PSNR & Minimum PSNR & Mean boundary RMSE \\ \midrule Mean geometry & $32.5922$ dB & $31.5426$ dB & $1.0261\times10^{-1}$ \\ Row/adjoint edge & $38.0604$ dB & $32.5913$ dB & $3.8393\times10^{-2}$ \\ Multi-ray FRI variable projection & $64.0770$ dB & $61.7645$ dB & $1.9800\times10^{-3}$ \\ True-geometry oracle & $97.9883$ dB & $83.1976$ dB & $0.0000$ \\ \bottomrule \end{tabular} \caption{Robustness sweep for multi-ray FRI edge refinement over $18$ cases: two boundary variants, three projection-angle budgets, and three noise levels. The multi-ray FRI geometry remains consistently high quality and closes most of the gap between crude adjoint peaks and the oracle.} \label{tab:haouchat-multiray-fri-robustness} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/haouchat_multiray_fri_robustness_outputs/haouchat_multiray_fri_robustness_sweep.png} \caption{Representative robustness-sweep preview for the base boundary at the largest projection budget.} \label{fig:haouchat-multiray-fri-robustness} \end{figure} \section{Sparse Stochastic Innovation Models} The sparse stochastic framework defines a process by \cite{unser2014sparse1,unser2014sparse2} \[ L\{s\}=w, \] where $L$ is a whitening operator and $w$ is white innovation noise. Gaussian $w$ produces dense least-sparse processes; non-Gaussian Levy noise produces sparse or impulsive innovations. The operator controls correlation and physics, while the Levy measure controls sparsity. The discrete-domain theory shows that matched B-spline filters convert continuous innovations into discrete generalized increments \cite{unser2014sparse2}. MAP and MMSE estimators for these priors are developed in \cite{bostan2013sparse,amini2013bayesian,kamilov2013mmse}. For OSNR, this means sparse parameters should not be arbitrary dense neural weights. They should be coefficient-domain innovations induced by the correct operator. \subsection{Controlled SPDE validation: advection--diffusion with Levy innovations} To convert the sparse stochastic theory into a PINN/weather-facing experiment, we implemented \texttt{apps\_industrial\_breakthrough/spde\_operator\_spline\_benchmark.py}. The controlled PDE is a periodic one-dimensional advection--diffusion--reaction model over a space--time block, \begin{equation} \mathcal{L}u = \left(\partial_t + a\partial_x-\nu\partial_{xx}+\lambda\right)u = w(x,t), \label{eq:spde-advection-diffusion} \end{equation} where $w$ is not restricted to be Gaussian. Following the Unser--Tafti sparse process model, the operator $\mathcal{L}$ fixes the correlation and propagation physics, while the innovation law determines the forcing morphology. We test four innovation profiles: a smooth periodic source, a Gaussian stochastic source, a compound-Poisson sparse impulse source, and a mixed weather-like source containing smooth waves, Gaussian background, and sparse jump events. On the periodic grid, \eqref{eq:spde-advection-diffusion} has the Fourier-domain symbol \begin{equation} \widehat{\mathcal{L}}(\omega_t,\omega_x) = j\omega_t + ja\omega_x+\nu\omega_x^2+\lambda, \end{equation} so the OSNR state-free solve is the diagonal complex division \begin{equation} \widehat{u}(\omega_t,\omega_x) = \frac{\overline{\widehat{\mathcal{L}}(\omega_t,\omega_x)}}{|\widehat{\mathcal{L}}(\omega_t,\omega_x)|^2+\epsilon} \widehat{w}(\omega_t,\omega_x). \label{eq:spde-fft-solve} \end{equation} This is the SPDE analogue of an operator-matched exponential spline solve: the Green structure is built into the inverse operator, and no neural coordinate residual or automatic-differentiation tape is required. A low-pass Fourier reconstruction is included as a spectral-bias baseline; it mimics what happens when a smooth model family cannot carry non-Gaussian sparse innovations. \begin{table}[h] \centering \scriptsize \begin{tabular}{lcccccc} \toprule Profile & OSNR PSNR & Low-pass PSNR & Innovation RMSE & Event error & Hard zeros & Solve time \\ \midrule Smooth periodic, $128^2$ & $138.0840$ dB & $85.4051$ dB & $1.2801\times10^{-4}$ & n/a & $0.00\%$ & $0.1383$ ms \\ Gaussian SPDE, $128^2$ & $90.9460$ dB & $38.7232$ dB & $2.9458\times10^{-4}$ & n/a & $0.00\%$ & $0.1370$ ms \\ Poisson sparse, $128^2$ & $69.9296$ dB & $34.4731$ dB & $3.6954\times10^{-3}$ & $0.0000$ px & $99.78\%$ & $0.1319$ ms \\ Mixed Levy weather, $128^2$ & $99.5862$ dB & $60.4309$ dB & $3.7530\times10^{-3}$ & $7.0438$ px & $0.00\%$ & $0.1414$ ms \\ Poisson sparse, low diffusion & $62.1205$ dB & $30.7575$ dB & $9.2609\times10^{-3}$ & $0.0000$ px & $99.41\%$ & $0.1353$ ms \\ Mixed Levy, low diffusion & $80.6083$ dB & $47.5993$ dB & $9.2813\times10^{-3}$ & $3.2817$ px & $0.00\%$ & $0.1472$ ms \\ Poisson sparse, $192^2$ & $71.4692$ dB & $34.3896$ dB & $2.4829\times10^{-3}$ & $0.0000$ px & $99.83\%$ & $0.5018$ ms \\ Mixed Levy weather, $192^2$ & $104.2661$ dB & $63.6906$ dB & $2.5052\times10^{-3}$ & $5.1575$ px & $0.00\%$ & $0.5070$ ms \\ \bottomrule \end{tabular} \caption{Controlled SPDE operator-spline benchmark for advection--diffusion with Gaussian and sparse Levy innovations. The OSNR solve is the direct FFT inversion of \eqref{eq:spde-fft-solve}; the low-pass row is a smooth spectral-bias baseline.} \label{tab:spde-operator-spline} \end{table} \begin{figure}[h] \centering \includegraphics[width=\linewidth]{../apps_industrial_breakthrough/spde_operator_spline_outputs_main/spde_operator_spline_profiles.png} \caption{SPDE profile comparison on the $128^2$ benchmark. Each row shows the target field, the OSNR state-free FFT reconstruction, and the low-pass smooth baseline. The sparse and mixed rows expose why a Gaussian/smooth-only surrogate is not enough for weather-like fronts and impulses.} \label{fig:spde-operator-spline} \end{figure} Table~\ref{tab:spde-operator-spline} gives the current interpretation. The state-free operator solve is essentially exact for all four innovation laws and remains below one millisecond even at $192^2$. The pure compound-Poisson case recovers event coordinates exactly at the tested grid resolutions, validating the sparse innovation view. The mixed case is more realistic and more difficult: the field reconstruction remains excellent, but raw top-$K$ event localization degrades because the smooth and Gaussian components overlap the sparse impulses in the recovered innovation. This is not a failure of the operator inverse; it identifies the next algorithmic requirement. A weather-grade OSNR solver should add the same sparse-plus-smooth oblique innovation sieve used elsewhere in this paper, but now applied to $\mathcal{L}u$ rather than to the field $u$ itself. We therefore added an explicit innovation-domain sieve to the same benchmark. Given the recovered innovation $\tilde{w}=\mathcal{L}\tilde{u}$, the sieve first estimates a smooth background $w_{\mathrm{sm}}=G_\sigma\ast \tilde{w}$ and then extracts sparse events from the residual \begin{equation} w_{\mathrm{sp}}=\mathcal{T}(\tilde{w}-w_{\mathrm{sm}}), \end{equation} where $\mathcal{T}$ is either a known-cardinality top-$K$ selector or an adaptive median-absolute-deviation threshold. This is not a field smoother; it acts after applying the physical operator and is therefore an innovation prior in the sense of Unser and Tafti. On the $128^2$ mixed Levy/weather case, raw top-$K$ localization has mean event error $7.0438$ px. The known-cardinality sieve reduces this to $0.0000$ px. The same result holds for the low-diffusion stress case, where raw localization is $3.2817$ px, and for the $192^2$ case, where raw localization is $5.1575$ px. The adaptive MAD sieve with threshold $8$ also recovers the mixed-weather events exactly without being told the number of events; it selects $36$ sparse sites in the mixed case and yields $0.0000$ px event error. In the pure Poisson case it selects a larger sparse support ($602$ sites at $128^2$) because the Gaussian smoothing residual leaves a local halo around each impulse, but nearest-event localization is still exact. Thus the next refinement is amplitude/support debiasing, not event detection. The practical weather implication is encouraging: OSNR can solve the stochastic PDE block globally and then separate sparse front/impulse innovations from smooth meteorological background in the physically meaningful residual domain. We also tested the immediate nonlinear extension in \texttt{apps\_industrial\_breakthrough/forced\_burgers\_spde\_benchmark.py}. The model is a periodically forced viscous Burgers equation, \begin{equation} u_t + uu_x-\nu u_{xx}=f_{\mathrm{smooth}}(x,t)+f_{\mathrm{sp}}(x,t), \end{equation} where $f_{\mathrm{sp}}$ is a sparse set of localized Gaussian events. A high-resolution spectral RK4 rollout is treated as the reference trajectory; compressed OSNR rollouts retain only a fixed number of Fourier/operator modes. The key inverse-problem distinction is that the forcing innovation must be estimated by applying the nonlinear physical operator to the observed trajectory, \begin{equation} \tilde f(x,t)=u_t+u u_x-\nu u_{xx}, \end{equation} not by thresholding the difference between a coarse rollout and the reference. The latter is mostly a truncation and phase-defect diagnostic. The former is the nonlinear analogue of the operator-domain innovation extraction used in the linear SPDE experiment. We therefore report both sparse support recovery and compact event-atom recovery. The point sieve thresholds $\tilde f-G_\sigma\ast \tilde f$ and measures whether each true event overlaps the recovered sparse support. The weak-form atom score integrates $\tilde f$ against anisotropic Gaussian test functions matched to the injected event scale and then applies non-maximum suppression. We also apply a local centroid debiasing step around each detected atom. This approximates \begin{equation} \eta_m=\langle \tilde f,\varphi_m\rangle, \end{equation} where $\varphi_m$ is a compact adjoint/test atom. On clean synthetic forcing, direct operator-domain support recovery is sharper than the weak score; the weak form is expected to become more useful once observations are noisy or irregular. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Coarse modes & Rollout PSNR & Support error & Atom error & Centroid error & Defect ratio \\ \midrule $8$ & $39.5805$ dB & $0.0000$ px & $1.0357$ px & $0.7917$ px & $0.2290$ \\ $18$ & $68.7114$ dB & $0.0000$ px & $1.0357$ px & $0.7917$ px & $0.0098$ \\ $32$ & $108.0643$ dB & $0.0000$ px & $1.0357$ px & $0.7917$ px & $0.0001$ \\ \bottomrule \end{tabular} \caption{Forced Burgers SPDE diagnostic after correcting the inverse-problem residual. Applying the nonlinear operator to the observed trajectory recovers every sparse forcing support location at the tested grid resolution. The atom-center error is about one pixel because the injected events are finite-width Gaussian blobs and overlapping events shift local maxima; local centroid debiasing reduces this to $0.7917$ px. The defect ratio reports the norm of the coarse-rollout phase/truncation defect relative to the physical innovation norm.} \label{tab:forced-burgers-spde} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/forced_burgers_spde_weak_outputs_main/forced_burgers_spde.png} \caption{Forced Burgers SPDE diagnostic at $18$ retained modes. The corrected operator innovation $\tilde f=u_t+u u_x-\nu u_{xx}$ exposes the sparse forcing structure directly. The support map recovers the event locations, while the weak-form score produces compact event atoms within about one pixel.} \label{fig:forced-burgers-spde} \end{figure} The conclusion is important for the weather/PINN program. Linear operator-matched SPDEs are already a home-turf win for OSNR: exact global solves, sparse Levy innovations, and sub-millisecond runtime. The corrected nonlinear Burgers diagnostic shows that sparse forcing can also be recovered when the physical operator is applied in the right domain. A denser stress case with $56$ injected events still yields $0.0000$ px support error and $0.8912$ px centroid error. The remaining bottleneck is not event detection but support and amplitude debiasing for finite-width/overlapping events, especially under noisy or partially observed fields. The next layer should estimate sparse innovations through an adjoint weak form, \begin{equation} \langle f,\varphi_m\rangle = \langle u_t+uu_x-\nu u_{xx},\varphi_m\rangle, \end{equation} with test functions $\varphi_m$ matched to the operator and the expected front scale, plus a local centroid/amplitude debiasing step. An operator-splitting scheme that alternates deterministic nonlinear advection with a sparse forcing inverse problem is the natural production path before claiming weather-grade nonlinear SPDE recovery. \paragraph{Direct PINN home-turf challenger.} We added a more direct PINN-facing control in \texttt{apps\_industrial\_breakthrough/pinn\_operator\_home\_turf\_challenger.py}. The benchmark is a periodic two-dimensional Helmholtz/Poisson problem, \begin{equation} (-\Delta+\lambda)u(x,y)=f(x,y), \end{equation} where $u$ is a mixed-frequency smooth field and $f$ is obtained by applying the known operator. The OSNR path solves the field by a single FFT-domain division. The baseline is a SIREN-style coordinate PINN trained with Adam on data samples and automatic-differentiation residual collocation. On the quality-first $128^2$ run with $\lambda=6$, the full OSNR solve reaches $141.7829$ dB PSNR and RMSE $1.9704\times10^{-7}$ in $0.1695$ ms on CPU. A compressed low-mode OSNR profile retaining only $8.3557\%$ of Fourier bins still reaches $136.3778$ dB in $0.2935$ ms. The SIREN PINN baseline, after $1800$ epochs, reaches only $21.1458$ dB and RMSE $2.1203\times10^{-1}$ after $109.112$ s. This is not a noisy external-data claim; it is a clean operator-known PINN control. It demonstrates the central home-turf point: when the differential operator and boundary topology are known, structural inversion gives both higher accuracy and roughly $6.44\times10^5$ lower training latency than residual-learning the same field. \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Profile & PSNR & RMSE & Time & Active coefficients \\ \midrule OSNR full spectral solve & $141.7829$ dB & $1.9704\times10^{-7}$ & $0.1695$ ms & $100.00\%$ \\ OSNR low-mode solve & $136.3778$ dB & $3.6712\times10^{-7}$ & $0.2935$ ms & $8.3557\%$ \\ SIREN PINN, $1800$ epochs & $21.1458$ dB & $2.1203\times10^{-1}$ & $109.112$ s & dense MLP \\ \bottomrule \end{tabular} \caption{Direct PINN home-turf challenger on a periodic $128^2$ Helmholtz/Poisson field. The OSNR rows are measured FFT/operator inversions; the SIREN PINN row is measured Adam training with autograd residual collocation.} \label{tab:pinn-home-turf} \end{table} \begin{figure}[h] \centering \includegraphics[width=\textwidth]{../apps_industrial_breakthrough/pinn_operator_home_turf_outputs_n128_e1800/pinn_operator_home_turf_panel.png} \caption{PINN home-turf visual panel. The full and low-mode OSNR inversions are visually indistinguishable from the target at the displayed scale, while the trained SIREN PINN remains visibly over-smoothed after the measured optimization budget.} \label{fig:pinn-home-turf} \end{figure} The scale follow-up at $256^2$ confirms that the coefficient fraction improves with resolution when the operator spectrum is compact. With the same low-mode budget, OSNR reaches $143.7870$ dB in $0.6588$ ms for the full solve, and $136.8069$ dB in $1.1293$ ms while retaining only $2.0889\%$ of Fourier bins. A $900$-epoch PINN baseline on the same field reaches $19.4427$ dB after $59.916$ s. At $512^2$, the same low-mode budget retains only $0.5222\%$ of Fourier bins and still reaches $136.4969$ dB in $3.5055$ ms; the full solve reaches $143.3450$ dB in $2.3255$ ms, while a $300$-epoch PINN baseline reaches $19.0573$ dB after $19.664$ s. The purpose of these rows is not to claim a universal neural-operator benchmark victory; they isolate the regime where PINN residual learning is structurally the wrong computational tool. We ran an additional observation-noise stress test to separate robust atom detection from brittle support thresholding. Gaussian observation noise is added to the trajectory before evaluating the nonlinear operator. At $0.1\%$ relative observation noise, raw pointwise support thresholding misses many events ($7.8618$ px support error), but ranked atom selection from the same operator residual remains accurate ($0.7801$ px after centroid refinement). Mild pre-operator smoothing restores support overlap ($0.0357$ px) but blurs atom centers ($1.4350$ px). At $0.5\%$ noise, the best tested atom setting uses $\sigma=0.75$ pre-smoothing and reaches $0.7975$ px centroid error, while binary support thresholding is unreliable. This confirms the correct noisy-weather design: detect a ranked set of operator-domain event atoms first, then run local amplitude/support debiasing rather than relying on a global hard threshold. \subsection{Operator-symbol identification by variable projection} The preceding SPDE experiments assume that the differential operator is known. The next weather/PINN question is whether OSNR can also learn a compact operator from data without falling back to a dense coordinate network. We therefore added \texttt{apps\_industrial\_breakthrough/operator\_pole\_identification\_benchmark.py}. The controlled model is the same advection--diffusion--reaction family \begin{equation} \left(\partial_t+a\partial_x-\nu\partial_{xx}+\lambda\right)u=w, \end{equation} but now the coefficients $(a,\nu,\lambda)$ are treated as unknown operator parameters. In Fourier space, \begin{equation} \widehat{w} - j\omega_t\widehat{u} = \left(ja\omega_x+\nu\omega_x^2+\lambda\right)\widehat{u}, \end{equation} so the unknown operator coefficients enter linearly once the observed field and innovation are transformed. OSNR therefore identifies the operator by one complex ridge least-squares solve over selected frequency bins, \begin{equation} \widehat{\theta} = \arg\min_{\theta=(a,\nu,\lambda)} \left\| \mathbf{D}(\widehat{u})\theta - \left(\widehat{w}-j\omega_t\widehat{u}\right) \right\|_2^2 +\epsilon\|\theta\|_2^2. \end{equation} This is a variable-projection step: linear field coefficients remain solved by the operator inverse, while the low-dimensional operator symbol is recovered directly from the data. A backpropagation baseline optimizes the same three parameters by Adam through the spectral residual for $800$ steps. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Profile & Method & $\hat a$ & $\hat\nu$ & $\hat\lambda$ & PSNR & Time \\ \midrule Clean, all bins & OSNR LS & $0.730001$ & $0.021008$ & $0.168636$ & $76.9176$ dB & $1.4462$ ms \\ Clean, band $24$ & OSNR LS & $0.729999$ & $0.021000$ & $0.170013$ & $116.0950$ dB & $0.2236$ ms \\ Clean & Adam residual & $0.729998$ & $0.037428$ & $0.038625$ & $17.9194$ dB & $206.4690$ ms \\ $0.5\%$ noise, band $24$ & OSNR LS & $0.729844$ & $0.019821$ & $0.366630$ & $40.2696$ dB & $0.1575$ ms \\ $0.5\%$ noise, band $12$ & OSNR LS & $0.729962$ & $0.020995$ & $0.171116$ & $56.3422$ dB & $0.1897$ ms \\ $0.5\%$ noise & Adam residual & $0.723947$ & $0.030782$ & $0.044012$ & $20.8344$ dB & $205.5888$ ms \\ \bottomrule \end{tabular} \caption{Operator-symbol identification for the advection--diffusion--reaction family with true parameters $(a,\nu,\lambda)=(0.73,0.021,0.17)$. The OSNR row uses a single complex least-squares solve in the Fourier/operator domain; the baseline uses iterative backpropagation through the same residual. Conservative spectral fitting bands suppress derivative-amplified observation noise.} \label{tab:operator-pole-identification} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/operator_pole_identification_outputs_noise005_band12/operator_pole_identification.png} \caption{Noisy operator-identification result at $0.5\%$ observation noise with a conservative fitting band. The recovered operator reconstructs the state at $56.3422$ dB after one structured coefficient solve.} \label{fig:operator-pole-identification} \end{figure} Table~\ref{tab:operator-pole-identification} is the first explicit operator-learning result. In the clean case, a frequency band of $24$ modes recovers all three coefficients to near machine precision and improves the reconstruction from $76.9176$ dB to $116.0950$ dB by avoiding ill-conditioned bins. With $0.5\%$ observation noise, fitting too many frequencies corrupts the reaction estimate because derivative operators amplify high-frequency noise. Tightening the band to $12$ modes restores the coefficients to sub-percent relative error and yields $56.3422$ dB, while the Adam residual baseline remains near $20.8$ dB after $800$ gradient steps. The lesson is directly relevant to weather data: unknown physics should be learned as a compact, stability-constrained operator symbol with explicit spectral/noise control, not as an unconstrained dense coordinate network. We then tested the harder field-only variant in \texttt{apps\_industrial\_breakthrough/blind\_operator\_sparsity\_identification.py}. Here $w$ is hidden: the search chooses the operator whose residual $\mathcal{L}_\theta u$ is most compressible as a low-pass smooth field plus a fixed number of sparse atoms. This is closer to unsupervised weather-model discovery, but it exposes an identifiability boundary. On four independent trajectories sharing the same true operator, the blind compressibility score selects $(\hat a,\hat\nu,\hat\lambda)=(0.91,0.021,0.26)$ instead of $(0.73,0.021,0.17)$, even though the sparse event locations are recovered exactly. A support-projected oracle that masks the true sparse event neighborhoods but does not know their amplitudes also fails to recover the reaction coefficient. The reason is structural: from $u$ alone, a wrong operator can be absorbed into a different smooth forcing background, so sparse-plus-smooth compressibility is not a unique operator identifier. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Setting & $\hat a$ & $\hat\nu$ & $\hat\lambda$ & Event error & Time \\ \midrule Blind sieve, $4$ clean trajectories & $0.910000$ & $0.021000$ & $0.260000$ & $0.0000$ px & $2596.29$ ms \\ Blind sieve, $4$ trajectories, $0.2\%$ noise & $0.910000$ & $0.009000$ & $0.260000$ & $0.0000$ px & $2594.86$ ms \\ Support-projected oracle, clean & $0.570370$ & $0.014276$ & $-14.202470$ & oracle support & $5.86$ ms \\ Support-projected oracle, $0.2\%$ noise & $0.314185$ & $0.000886$ & $1.900169$ & oracle support & $5.96$ ms \\ \bottomrule \end{tabular} \caption{Blind operator-discovery diagnostic. Sparse event geometry can be recovered from $\mathcal{L}_\theta u$, but field-only sparse-plus-smooth compressibility does not uniquely identify the true operator because operator mismatch can be reinterpreted as smooth forcing.} \label{tab:blind-operator-sparsity} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/blind_operator_sparsity_outputs_support_projected/blind_operator_sparsity_identification.png} \caption{Blind operator-sparsity diagnostic. The selected residual preserves sparse event locations but corresponds to the wrong operator, demonstrating that fully blind field-only operator discovery needs additional physical anchors.} \label{fig:blind-operator-sparsity} \end{figure} This negative result is useful. It says the breakthrough lane is not arbitrary unsupervised PDE discovery from a single scalar field. The credible path is semi-blind operator learning: use measured innovations, multiple observed state channels, conservation laws, boundary/flux constraints, or assimilation windows to anchor the smooth forcing ambiguity, then recover the compact operator symbol by structured least squares or variable projection. The first semi-blind anchor test follows this prescription. In \texttt{apps\_industrial\_breakthrough/anchored\_operator\_identification.py}, only a random subset of the forcing samples is revealed. The field $u$ is observed everywhere, but the operator coefficients are fitted from the pointwise equations \begin{equation} u_t(t_i,x_i)+a u_x(t_i,x_i)-\nu u_{xx}(t_i,x_i)+\lambda u(t_i,x_i)=w(t_i,x_i) \end{equation} at the anchor sites. All derivatives are evaluated analytically by spectral/operator columns, and the three unknown coefficients are recovered by a real ridge least-squares solve. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Setting & Anchors & $\hat a$ & $\hat\nu$ & $\hat\lambda$ & PSNR \\ \midrule Clean, $0.25\%$ anchors & $41$ & $0.729829$ & $0.021105$ & $0.143701$ & $49.6630$ dB \\ Clean, $0.50\%$ anchors & $82$ & $0.729892$ & $0.021072$ & $0.159291$ & $58.4940$ dB \\ Clean, $5.00\%$ anchors & $819$ & $0.729956$ & $0.021002$ & $0.172572$ & $69.6683$ dB \\ $0.2\%$ noise, no denoise, $5.00\%$ anchors & $819$ & $0.732336$ & $0.008352$ & $2.355422$ & $29.6130$ dB \\ $0.2\%$ noise, band $12$, $5.00\%$ anchors & $819$ & $0.729782$ & $0.020778$ & $0.183559$ & $53.7874$ dB \\ \bottomrule \end{tabular} \caption{Semi-blind operator identification from sparse forcing anchors. Clean operator recovery is accurate with very few forcing samples. Under observation noise, derivative columns require spectral denoising; with a band-$12$ field prefilter, $5\%$ anchors recover the operator to below $8\%$ worst relative error and reconstruct the state at $53.7874$ dB.} \label{tab:anchored-operator-identification} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/anchored_operator_identification_outputs_noise002_denoise12_moreanchors/anchored_operator_identification.png} \caption{Semi-blind noisy operator identification with sparse forcing anchors. A small set of pointwise forcing measurements breaks the field-only ambiguity exposed in Table~\ref{tab:blind-operator-sparsity}.} \label{fig:anchored-operator-identification} \end{figure} Table~\ref{tab:anchored-operator-identification} is a more realistic weather/PINN direction than fully blind scalar discovery. It shows that a small number of physical anchors can make compact operator identification well-posed again. The nonmonotone noisy rows also identify the next engineering layer: anchors should be selected by leverage or derivative-energy criteria rather than uniformly at random. We tested the simplest version of that idea by selecting anchors with the largest normalized derivative-column energy. This naive leverage rule is not sufficient. In the clean case it reaches only $52.1897$ dB at $5\%$ anchors, below the random-anchor $69.6683$ dB result. With $0.2\%$ observation noise and band-$12$ denoising, it degrades to $37.4777$ dB at $5\%$ anchors because high-leverage points are also the points where derivative noise is most amplified. The active-anchor rule must therefore combine derivative leverage with noise sensitivity and spatial diversity; selecting the largest rows of the design matrix is too brittle. A follow-up robustification adds a trimmed ridge solve: after the first anchor fit, the largest pointwise residuals are discarded and the operator is refit on the lowest-residual fraction. Random anchors with trimming do not improve the best $5\%$ noisy row, but diverse leverage plus trimming uncovers a useful low-anchor operating point. With only $0.5\%$ anchors under $0.2\%$ observation noise, diverse-trimmed anchors estimate $(a,\nu,\lambda)=(0.731760,0.020692,0.175461)$, corresponding to only $3.2124\%$ worst relative parameter error and $49.5094$ dB reconstruction. The best state PSNR still comes from denser random anchors, but the active-trimmed result shows that carefully chosen anchors can reduce physical measurements by an order of magnitude while preserving an accurate compact operator. \subsection{Nonlinear Burgers operator identification and held-out forecasting} The next SOTA-facing PINN target is nonlinear forecasting rather than static reconstruction. We implemented \texttt{apps\_industrial\_breakthrough/burgers\_operator\_identification\_forecast.py}, which treats viscous Burgers dynamics \begin{equation} u_t + c\,u u_x = \nu u_{xx} \end{equation} as a compact operator-identification problem. From the observed training window, OSNR forms analytic derivative columns $(u_t,uu_x,u_{xx})$ and solves the two unknown coefficients $(c,\nu)$ by a tiny ridge system over sparse anchor samples. The learned operator is then rolled forward over the held-out future window. This is the nonlinear analogue of the semi-blind anchor experiments above, but the validation target is future prediction, which is the quantity that PINNs and neural operators usually report. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Setting & Anchors & $\hat c$ & $\hat\nu$ & Future PSNR & ID time \\ \midrule Clean, $0.25\%$ anchors & $35$ & $1.003308$ & $0.0045239$ & $70.7236$ dB & $0.4064$ ms \\ Clean, $0.50\%$ anchors & $69$ & $0.998324$ & $0.0044914$ & $77.6890$ dB & $0.2062$ ms \\ Clean, $1.00\%$ anchors & $138$ & $0.999828$ & $0.0044959$ & $88.6117$ dB & $0.1814$ ms \\ Adam residual, clean & all & $0.997542$ & $0.0044901$ & $74.8301$ dB & $82.15$ ms \\ $0.1\%$ noise, $0.50\%$ anchors & $69$ & $0.995605$ & $0.0044539$ & $66.3736$ dB & $0.2458$ ms \\ $0.2\%$ noise, $1.00\%$ anchors & $138$ & $0.995355$ & $0.0044588$ & $66.8131$ dB & $0.3160$ ms \\ Adam residual, $0.2\%$ noise & all & $0.979668$ & $0.0044440$ & $56.8246$ dB & $81.62$ ms \\ Wrong prior & n/a & $0.75$ & $0.00225$ & $30.7955$ dB & n/a \\ \bottomrule \end{tabular} \caption{Nonlinear Burgers operator identification and held-out forecasting. True parameters are $(c,\nu)=(1.0,0.0045)$, the training window is $45\%$ of the timeline, and the future window contains the remaining $90$ frames. OSNR identifies the operator in sub-millisecond time from sparse anchors and forecasts the future without a neural training loop.} \label{tab:burgers-operator-id-forecast} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/burgers_operator_identification_outputs_noise002_modes28_sub16/burgers_operator_identification_forecast.png} \caption{Noisy Burgers operator-ID forecast at $0.2\%$ observation noise. OSNR recovers the nonlinear operator from $1\%$ training-window anchors and forecasts the held-out future at $66.8131$ dB.} \label{fig:burgers-operator-id-forecast} \end{figure} This is the strongest nonlinear PINN-facing result so far. The clean $1\%$ anchor row reaches $88.6117$ dB future PSNR with a $0.1814$ ms identification solve, while the Adam residual fit is about $450\times$ slower and reaches only $74.8301$ dB. Under $0.2\%$ observation noise, OSNR still reaches $66.8131$ dB from $1\%$ anchors, outperforming the Adam residual fit by about $10$ dB. The wrong-prior row shows that the forecast is not trivially easy: incorrect physics collapses to $30.7955$ dB. This is the first result that directly combines nonlinear coefficient discovery, held-out forecasting, sparse physical measurements, and a clear optimization-speed gap. \subsection{Coupled nonlinear shallow-water operator forecasting} The scalar Burgers forecast is a necessary control, but the weather/PINN claim requires a coupled multi-field nonlinear system. We therefore implemented \texttt{apps\_industrial\_breakthrough/nonlinear\_shallow\_water\_operator\_forecast.py}. The state is $\mathbf q=(\eta,u,v)$ and the operator family is \begin{align} \eta_t &= -H(u_x+v_y)-\beta\,\nabla\cdot(\eta(u,v))+\mu_h\Delta\eta,\\ u_t &= -g\eta_x+f v-r u-\beta(u u_x+v u_y)+\nu\Delta u,\\ v_t &= -g\eta_y-f u-r v-\beta(u v_x+v v_y)+\nu\Delta v. \end{align} The unknown physical vector is \[ \theta=(H,\beta,g,f,r,\nu,\mu_h), \] covering mean depth, nonlinear transport strength, gravity, Coriolis coupling, damping, momentum viscosity, and height diffusion. OSNR forms the full analytic derivative library from the observed training window and solves a scaled sparse-anchor linear system for all seven coefficients at once. The fitted operator is then advanced over the held-out future window using the same spectral RK4 physics core. The comparison baseline fits the same residual equations by Adam over all training rows, so the quality comparison is not against a weak interpolant but against the standard differentiable residual-minimization path used by PINN-style methods. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Setting & Anchors & $\hat\beta$ & $\hat g$ & Max param. err. & Future PSNR & ID time \\ \midrule $48^2$, clean, $0.25\%$ & $795$ & $0.7178$ & $0.8590$ & $1.8214\%$ & $68.7348$ dB & $14.24$ ms \\ $48^2$, clean, $1.00\%$ & $3180$ & $0.7181$ & $0.8590$ & $1.6975\%$ & $68.7314$ dB & $13.24$ ms \\ Adam residual, $48^2$ clean & all & $0.7186$ & $0.8590$ & $1.8572\%$ & $68.6608$ dB & $830.30$ ms \\ $48^2$, $0.1\%$ noise & $795$ & $0.7091$ & $0.8594$ & $21.347\%$ & $66.3867$ dB & $16.93$ ms \\ Adam residual, $0.1\%$ noise & all & $0.7171$ & $0.8587$ & $17.495\%$ & $66.0641$ dB & $1093.02$ ms \\ $48^2$, $0.2\%$ noise & $795$ & $0.7218$ & $0.8599$ & $34.417\%$ & $65.6597$ dB & $17.96$ ms \\ Adam residual, $0.2\%$ noise & all & $0.7156$ & $0.8585$ & $32.218\%$ & $63.4251$ dB & $1095.57$ ms \\ $64^2$, clean, $0.10\%$ & $713$ & $0.7191$ & $0.8594$ & $1.8368\%$ & $72.6945$ dB & $24.00$ ms \\ Adam residual, $64^2$ clean & all & $0.7191$ & $0.8593$ & $1.1991\%$ & $72.5755$ dB & $1103.88$ ms \\ Wrong prior, $64^2$ & n/a & $0.3240$ & $1.0750$ & n/a & $27.7993$ dB & n/a \\ \bottomrule \end{tabular} \caption{Coupled nonlinear shallow-water operator identification and held-out forecasting. True parameters are $(H,\beta,g,f,r,\nu,\mu_h)=(1.0,0.72,0.86,0.58,0.065,0.006,0.004)$. OSNR uses sparse derivative anchors; Adam optimizes the same residual over all training rows.} \label{tab:nonlinear-shallow-water-forecast} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/nonlinear_shallow_water_operator_forecast_outputs_121x64_clean/nonlinear_shallow_water_operator_forecast.png} \caption{Coupled nonlinear shallow-water forecast at $121\times64^2$. With only $0.1\%$ sparse derivative anchors, OSNR identifies the seven-parameter nonlinear operator and forecasts the held-out future at $72.6945$ dB. The wrong-prior forecast falls to $27.7993$ dB, confirming that the high score is not a trivial smoothness artifact.} \label{fig:nonlinear-shallow-water-forecast} \end{figure} This is the first multi-field nonlinear weather-core forecast result in the project. It preserves the key advantage seen in Burgers: the residual landscape can be collapsed into a small structured operator solve instead of optimized by thousands of neural/PINN gradient steps. On the $64^2$ run, OSNR uses only $713$ anchor equations out of the full derivative library and identifies the operator in $24.00$ ms, while the Adam residual fit takes $1103.88$ ms. Both methods converge to similar coefficients in the clean case because the library is correct, but OSNR reaches the solution in one scaled linear solve with about a $46\times$ identification-speed advantage and no neural training loop. Under observation noise, derivative bias still affects the weak damping/diffusion terms, but the future forecast remains above $65$ dB and stays ahead of Adam in the tested $0.2\%$ setting. The next moonshot is therefore not another scalar PDE; it is sparse/partial observation data assimilation for this same coupled nonlinear operator family. \subsection{Adaptive sparse sensors for nonlinear shallow-water assimilation} We then cross-pollinated the weather station-placement result with the coupled nonlinear shallow-water core. The runner \texttt{apps\_industrial\_breakthrough/nonlinear\_shallow\_water\_adaptive\_sensor\_assimilation.py} keeps the same no-backprop pipeline: sparse sensors reconstruct the observed training window by a closed-form Fourier-dictionary ridge solve, the seven-parameter nonlinear operator is identified by sparse least squares, and the state is rolled into the held-out future. The only changed variable is where the sparse sensors are placed. Sensor policies are computed from the training window only, never from held-out future frames. We compare random points, raw training-window variance/gradient/leverage scores, lattice-plus-score hybrids, residual-innovation hybrids, and centered/phase-shifted coverage policies. \begin{table}[H] \centering \small \begin{tabular}{lcccc} \toprule Observation setting & Modes & Sensor policy & Sensors & Future PSNR \\ \midrule Clean & $4$ & lattice $3\%$ & $123$ & $19.842$ dB \\ Clean & $4$ & lattice+hybrid $10\%$ budget, $3\%$ total & $123$ & $19.822$ dB \\ Clean & $3$ & centered lattice $1.5\%$ & $61$ & $19.852$ dB \\ Clean & $3$ & centered lattice $2\%$ & $82$ & $19.942$ dB \\ Clean & $3$ & offset-best lattice $2\%$ & $82$ & $19.958$ dB \\ Clean & $3$ & random $20\%$ & $819$ & $19.920$ dB \\ $1\%$ sensor noise & $3$ & centered lattice $2\%$ & $82$ & $19.948$ dB \\ $1\%$ sensor noise & $3$ & random $20\%$ & $819$ & $19.906$ dB \\ \bottomrule \end{tabular} \caption{Adaptive sparse-sensor placement for nonlinear shallow-water assimilation and forecasting. All rows use the same $121\times64^2$ trajectory, closed-form sparse-window assimilation, sparse least-squares operator identification, and no backpropagation. Mode $3$ centered/offset coverage reaches dense-random forecast quality with $10\times$ fewer observations. The offset-best policy chooses the best phase among $16$ centered lattice shifts by training-window reconstruction PSNR only.} \label{tab:nonlinear-shallow-water-adaptive-sensors} \end{table} \begin{figure}[H] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/nonlinear_shallow_water_adaptive_sensor_assimilation_outputs_full_modes3_offset_lattice/adaptive_sensor_panel.png} \caption{Adaptive shallow-water sparse-sensor forecast panel for the mode-$3$ coverage-geometry run. Centered/offset lattice rows preserve the large-scale future height field with $2$--$3\%$ sensors, while dense random placement needs about $20\%$ sensors to reach the same forecast band.} \label{fig:nonlinear-shallow-water-adaptive-sensors} \end{figure} This result is useful because it is positive and diagnostic. The naive high-information policies are not winners: variance, gradient, and hybrid-diverse placement overconcentrate sensors in active regions and can make the Fourier reconstruction ill-conditioned. The follow-up lattice-plus-information experiment confirmed the same boundary: at $3\%$ total sensors with mode $4$, lattice+gradient, lattice+hybrid, and lattice+residual $10\%$ allocation reach $19.786$, $19.822$, and $19.807$ dB, all below the pure lattice row at $19.828$ dB. The real improvement is coverage geometry plus basis order. Mode $3$ is the sparse bias-variance sweet spot; modes $1$--$2$ underfit and modes $5$--$6$ are underconstrained at low sensor counts. A centered lattice at $2\%$ sensors reaches $19.942$ dB clean future PSNR, above the same-run $20\%$ random reference at $19.920$ dB; with $1\%$ sensor noise, the same $2\%$ centered lattice row reaches $19.948$ dB versus noisy random $20\%$ at $19.906$ dB. The low-count sweep shows the transition: $0.5\%$ centered sensors fail ($15.628$ dB), $1\%$ is not yet dense-random quality ($19.158$ dB), $1.5\%$ approaches it ($19.852$ dB), and $2\%$ crosses it. We then tested whether this was a single-trajectory phase artifact. The script now exposes initial roll and amplitude controls, and the offset-best policy chooses the best of $16$ lattice phases by training-window reconstruction only. Across four robustness variants, offset-best $2\%$ sensors remains at or above random $20\%$: roll $(7,11)$ gives $19.946$ versus $19.942$ dB, roll $(13,5)$ gives $19.955$ versus $19.919$ dB, amplitude scale $1.25$ gives $19.948$ versus $19.894$ dB, and a changed dynamics profile $(\beta,g,f)=(0.9,0.78,0.45)$ gives $19.981$ versus $19.932$ dB. Thus adaptive station placement is not only a terminal weather-assimilation trick. In a coupled nonlinear weather-core forecast, a coverage-aware sensor topology can reduce observations by $10\times$ while preserving future forecast quality, using only algebraic assimilation and operator identification. Finally, we removed the known-family assumption after sparse sensing. The runner \texttt{apps\_industrial\_breakthrough/nonlinear\_shallow\_water\_sparse\_sensor\_library\_discovery.py} first reconstructs the training window from sparse sensors, then fits the $24$-column shallow-water library by sequential thresholded least squares, and forecasts from the assimilated last state. This is a harder test because the solver must reject decoy columns and no longer receives the seven-parameter operator family. \begin{table}[H] \centering \scriptsize \setlength{\tabcolsep}{3pt} \begin{tabular}{@{}p{0.17\linewidth}p{0.13\linewidth}c c c c p{0.15\linewidth}@{}} \toprule Observation setting & Policy & Sensors & Threshold & Support $(TP,FP,FN)$ & Future PSNR & Reference \\ \midrule Clean & offset-best $2\%$ & $82$ & $0.003$ & $(8,0,5)$ & $19.960$ dB & known-family $19.962$ dB \\ Clean & random $20\%$ & $819$ & $0.005$ & $(7,0,6)$ & $19.888$ dB & known-family $19.902$ dB \\ $1\%$ sensor noise & offset-best $2\%$ & $82$ & $0.003$ & $(8,0,5)$ & $19.965$ dB & known-family $19.964$ dB \\ $1\%$ sensor noise & random $20\%$ & $819$ & $0.003$ & $(7,0,6)$ & $19.885$ dB & known-family $19.912$ dB \\ \bottomrule \end{tabular} \caption{Sparse-sensor governing-equation discovery after Fourier assimilation. The discovered support is counted against the $13$ true library columns and $11$ decoys. Offset-best $2\%$ sensors recover an $8$-term true subset with zero decoys and match or exceed dense-random $20\%$ forecast quality.} \label{tab:nonlinear-shallow-water-sparse-sensor-library} \end{table} The sparse-library result changes the interpretation. The $2\%$ offset-best row does not fully recover all weak nonlinear/damping terms, but it recovers the dominant conservative, pressure, Coriolis, and diffusion operators with no decoys and forecasts within about $0.002$ dB of the known-family coefficient fit. A lower threshold $0.001$ recovers $9$ true terms with one decoy and reaches $19.964$ dB clean, but the zero-decoy $8$-term threshold is the cleaner scientific claim. Thus the sensor topology is not merely helping a fixed PDE prior; it preserves enough operator information for sparse governing-equation discovery from partial observations. The next step was to remove a flaw in the pointwise library: it thresholds every equation-specific term and every decoy on the same normalized scale, even though the physically meaningful shallow-water operators are shared typed groups. We therefore added a weak-form typed library in \texttt{nonlinear\_shallow\_water\_sparse\_sensor\_weakform\_discovery.py}. Each window enforces \[ q(t_b)-q(t_a)\approx \int_{t_a}^{t_b}\mathcal{L}_j(q(t))\,dt \] and projects the balance onto low Fourier test modes. The $13$ true pointwise terms are then tied into seven shared physical groups $(H,\beta,\mu_h,g,f,r,\nu)$, while the $11$ nuisance columns receive a larger typed-selection threshold. This is not a neural loss or a backward pass: it is a weak-form operator balance followed by weighted sequential thresholded least squares. \begin{table}[H] \centering \scriptsize \setlength{\tabcolsep}{3pt} \begin{tabular}{@{}p{0.17\linewidth}p{0.14\linewidth}c c c c p{0.17\linewidth}@{}} \toprule Observation setting & Policy & Sensors & $(\tau,\lambda_{\rm decoy})$ & Support $(TP,FP,FN)$ & Future PSNR & Reference \\ \midrule Clean & random $2\%$ & $82$ & $(2{\times}10^{-6},50)$ & $(7,1,0)$ & unstable & known-family $16.055$ dB \\ Clean & offset-best $2\%$ & $82$ & $(5{\times}10^{-7},20)$ & $(7,0,0)$ & $19.961$ dB & weak dense $19.961$ dB \\ Clean & random $20\%$ & $819$ & $(10^{-6},50)$ & $(7,0,0)$ & $19.909$ dB & weak dense $19.909$ dB \\ $1\%$ sensor noise & offset-best $2\%$ & $82$ & $(5{\times}10^{-7},20)$ & $(7,0,0)$ & $19.956$ dB & weak dense $19.956$ dB \\ $1\%$ sensor noise & random $20\%$ & $819$ & $(10^{-6},50)$ & $(7,0,0)$ & $19.910$ dB & weak dense $19.910$ dB \\ \bottomrule \end{tabular} \caption{Typed weak-form sparse-sensor governing-equation discovery. Support is counted over seven shared physical operator groups and eleven decoys. The decoy multiplier $\lambda_{\rm decoy}$ applies only to nuisance columns. Offset-best $2\%$ sensors recover the complete seven-group shallow-water operator with zero decoys and match the dense weak-form solve, while using $10\times$ fewer observations than random $20\%$.} \label{tab:nonlinear-shallow-water-typed-weak-sensor-discovery} \end{table} This closes the support-recovery gap left by Table~\ref{tab:nonlinear-shallow-water-sparse-sensor-library}. With typed weak-form rows, offset-best $2\%$ sensors recover all seven physical groups with zero false positives in both the clean and $1\%$ sensor-noise settings. The same row is also forecast-competitive: $19.961$ dB clean and $19.956$ dB noisy, above the corresponding typed random-$20\%$ rows ($19.909$ and $19.910$ dB). Random $2\%$ still fails despite selecting most physical groups, which confirms that the result is not merely a threshold artifact. The station geometry must preserve a well-conditioned weak operator balance; once it does, typed OSNR selection can recover the full coupled shallow-water operator from sparse partial observations without backpropagation. The typed weak-form result also passed the first robustness sweep. At $1\%$ offset sensors, the selector already recovers $(7,0,0)$ support but only reaches $19.505$ dB, so full support and forecast-quality crossing are separate requirements. At $2\%$ offset sensors, all four trajectory variants recover $(7,0,0)$: roll $(7,11)$ gives $19.926$ dB, roll $(13,5)$ gives $19.948$ dB, initial scale $1.25$ gives $19.946$ dB, and changed dynamics $(\beta,g,f)=(0.9,0.78,0.45)$ gives $19.978$ dB. The random-$20\%$ typed weak-form rows for the same variants are $19.942$, $19.916$, $19.944$, and $19.977$ dB, respectively. Thus complete support recovery is robust in the tested variants; the $10\times$ forecast-quality advantage holds in three of four variants and narrowly fails on roll $(7,11)$. We then tested whether the same typed weak-form mechanism learns a reusable operator rather than a trajectory-specific correction. The multi-trajectory runner trains one shared operator from sparse-observed variants \{base, roll $(7,11)$, scale $1.25$\} and forecasts unseen roll variants $(13,5)$ and $(5,17)$ from their held-out states. Offset-best $2\%$ sensors recover full support and reach a mean unseen-trajectory PSNR of $51.100$ dB; random $2\%$ is unstable even with nearly full support. Dense random $20\%$ also recovers full support and reaches $52.232$ dB, while offset-best $20\%$ reaches $62.902$ dB. This is the first sparse-observation cross-trajectory operator-learning result in this section. It is not a $10\times$ dense-quality win at $2\%$ observations, but it shows that the typed OSNR weak-form solver can learn a shared coupled operator from partial observations and transfer it to unseen initial conditions without backpropagation. A follow-up sensor-fraction and placement sweep showed that the cross-trajectory coefficient bottleneck is primarily geometric. With offset-best placement and the same typed solver, $3\%$ sensors already recover $(7,0,0)$ and reach $53.282$ dB, exceeding the random-$20\%$ result; $5\%$ reaches $57.847$ dB; $10\%$ drops to $52.865$ dB; and $20\%$ reaches $62.902$ dB. Thus adding sensors is not monotone unless the station geometry remains well conditioned for the weak-form operator rows. A placement sweep at $3$--$10\%$ found that pointwise saliency policies (variance, gradient, and hybrid-diverse additions) consistently introduce decoy groups and degrade transfer. The best sparse clean-support result is a centered lattice with $5\%$ sensors and a stronger decoy multiplier: it recovers $(7,0,0)$ and reaches $58.829$ dB on the two unseen trajectories. If the two derivative-decoy groups are allowed, the same $5\%$ centered lattice reaches $60.775$ dB, but we treat this as a numerical correction rather than a clean governing-equation discovery. The practical conclusion is that operator-identifiability-balanced station geometry is more important than generic high-activity station placement. We then made the geometry test explicit by adding row-conditioning and fixed lattice-phase sweeps. Unit-design, unit-joint, and clipped row normalizations all made the solver worse: they activated most decoys and collapsed transfer to roughly $22$--$30$ dB. This negative result is important because the weak-form row magnitudes carry physical operator information; flattening them destroys the balance rather than improving conditioning. In contrast, fixed quarter-phase lattice placement is a productive control variable. At $5\%$ sensors on the three-training-trajectory protocol, the best clean fixed phase $(0.75,0.75)$ reaches $60.365$ dB with exact $(7,0,0)$ support, and the result validates on fresh roll and scale variants with mean $60.365$ dB. Expanding the training set to eight sparse-observed trajectories raises the same clean $5\%$ phase result to $60.945$ dB. Most importantly, a focused sensor curve with this phase shows that $7.5\%$ sensors per training trajectory recover exact support and reach $65.072$ dB on four fresh test variants, exceeding the same-protocol $20\%$ phase reference of $64.149$ dB. The curve remains nonmonotone: $10\%$ drops to $55.014$ dB and $15\%$ to $59.252$ dB. Thus the current lesson is sharper than ``more sensors'': sparse OSNR operator learning can beat denser observation budgets when station geometry is phase-balanced for the weak operator, but station-count increases can still harm coefficient estimation if they alias the weak-form rows. We also audited whether the phase can be selected without looking at the final test variants. Simple training-only proxies failed: physical-column condition number, physical--decoy coherence, dense weak residual, leave-one-training-trajectory weak residual, leave-one theta stability, and an inner training-window rollout score did not rank the best phases. A standard validation split, however, does. Selecting the phase on validation variants \{roll $(3,9)$, scale $0.75$\} chooses $(0.25,0.75)$, which then transfers to disjoint test variants \{roll $(11,4)$, scale $1.40$, roll $(19,2)$, scale $1.30$\}. The validation-selected $7.5\%$ phase recovers exact support and reaches $65.161$ dB on that disjoint test set, while the same phase with $20\%$ sensors reaches $64.283$ dB. This is the cleanest current sparse cross-trajectory result: the station geometry is selected on validation data, the test variants are unseen, and the learned seven-group operator still beats the denser observation budget without backpropagation. Finally, we tested whether the nonmonotone phase-lattice curve could be repaired by replacing the lattice with low-discrepancy or jittered station families. It could not. On the same validation-selected test protocol, Sobol phase stations activated $7$--$10$ decoys and produced unstable forecasts at $7.5\%$, $10\%$, and $15\%$ sensors. Jittered phase lattices were stable but much weaker: $46.774$ dB at $7.5\%$, $44.902$ dB at $10\%$, and $58.205$ dB at $15\%$. The unjittered phase lattice remains the best clean geometry, with $65.161$ dB at $7.5\%$. Thus the current station rule is not generic space filling; it is a Fourier-compatible phase-balanced sampling rule. We then promoted the validation split from phase selection to joint phase/count selection. The training set stayed fixed at eight sparse-observed variants \{base, roll $(7,11)$, scale $1.25$, roll $(13,5)$, roll $(5,17)$, roll $(2,19)$, scale $0.90$, scale $1.10$\}. The validation variants were again roll $(3,9)$ and scale $0.75$. The grid searched lattice phases $(0.25,0.75)$, $(0.75,0.25)$, $(0,0.5)$, $(0.75,0.75)$, $(0,0.75)$, and $(0.25,0.5)$ at sensor fractions $5\%$, $6.25\%$, $7.5\%$, $8.75\%$, $10\%$, $12.5\%$, and $15\%$, with threshold $5\times10^{-7}$, decoy multiplier $20$, no row normalization, window $7$, stride $2$, and a $75\%$ inner training-window diagnostic split. Validation selected the $8.75\%$ lattice with phase $(0.75,0.75)$: it recovered exact $(7,0,0)$ support and reached $67.996$ dB on the two validation variants. Without changing any hyperparameter, the selected row transferred to disjoint test variants \{roll $(11,4)$, scale $1.40$, roll $(19,2)$, scale $1.30$\}, reaching $67.841$ dB with exact support. Same-run references were $65.432$ dB at $6.25\%$, $65.161$ dB for the previous $7.5\%$ phase $(0.25,0.75)$ row, and $64.283$ dB for the same $20\%$ phase $(0.25,0.75)$ row. Thus held-out validation can now select both observation count and phase, and the selected sparse geometry uses only $358$ stations per training trajectory while outperforming $819$-station dense-phase references. A final fine phase/count refinement around this winner exposed a sharper resonance. We searched $8.125\%$, $8.4375\%$, $8.75\%$, $9.0625\%$, and $9.375\%$ sensors with phases $(0.62,0.62)$, $(0.62,0.75)$, $(0.75,0.62)$, $(0.75,0.75)$, $(0.75,0.87)$, $(0.87,0.75)$, $(0.87,0.87)$, $(0.62,0.87)$, and $(0.87,0.62)$. The adjacent count bands $8.125\%$ and $8.4375\%$ were poor despite exact support, reaching only about $53$--$54$ dB; $9.0625\%$ recovered to $68.520$ dB at phase $(0.75,0.75)$, but the validation winner was again $8.75\%$, now with phase $(0.75,0.87)$. This row reached $70.774$ dB on validation and transferred to the disjoint test variants at $69.803$ dB, with exact $(7,0,0)$ support and a $1.31$ ms sparse solve. The same fine-test run reproduced the old $8.75\%$ phase $(0.75,0.75)$ result at $67.841$ dB and showed that moving the winning phase to $9.0625\%$ drops to $65.941$ dB. The result is therefore not a generic phase preference. It is a count-specific Fourier sampling geometry that materially improves coefficient accuracy while keeping the observation budget at $358$ stations per training trajectory. To check whether the fine geometry was overfitting the two validation variants, we ran a broader fresh-variant audit with roll shifts $(1,23)$, $(23,1)$, $(31,17)$, $(17,31)$ and amplitude scales $0.60$, $1.60$, $0.50$, and $1.75$. The selected $8.75\%$ phase $(0.75,0.87)$ row reached $70.166$ dB across these eight variants with exact support. Same-count controls were $67.882$ dB for phase $(0.75,0.75)$ and $63.047$ dB for phase $(0.25,0.75)$, while $20\%$ references reached only $64.044$, $64.104$, and $65.107$ dB for the three tested phases. This robustness audit strengthens the interpretation: the selected sparse station geometry generalizes across unseen roll and amplitude perturbations and beats substantially denser station budgets because it better identifies the weak operator coefficients, not because it sees more observations. Because the resonance was phase-sharp, we then ran a local phase-only refinement at the fixed $8.75\%$ count. The validation grid swept $x$ phases $0.70,0.72,0.75,0.78,0.80$ and $y$ phases $0.84,0.87,0.90,0.93$ around the previous winner. Most rows were much weaker even with exact support; for example $x=0.80$ remained below $60$ dB and $(0.75,0.93)$ dropped to $66.390$ dB. The validation winner was $(0.75,0.90)$ at $73.568$ dB. Tested on the union of the four disjoint variants and the eight broad-audit variants, this row reached $73.379$ dB with exact support and a $1.17$ ms solve. On the same $12$-variant audit, $(0.75,0.87)$ reached $70.045$ dB and $(0.75,0.75)$ reached $67.868$ dB. This is the current best clean shallow-water result: validation-selected sparse station geometry with $358$ observations per training trajectory beats both same-count neighboring phases and all tested $819$-station references by a large margin. One more one-percent refinement around $(0.75,0.90)$ saturated rather than improved the result. Sweeping $x\in\{0.73,0.74,0.75,0.76,0.77\}$ and $y\in\{0.88,0.89,0.90,0.91,0.92\}$ at the same $8.75\%$ count again selected $(0.75,0.90)$; $(0.75,0.91)$ tied it because the rounded station set is effectively equivalent. Nearby rows drop quickly: $(0.75,0.89)$ gives $72.614$ dB, $(0.75,0.92)$ gives $68.431$ dB, $x=0.76$--$0.77$ with $y=0.90$--$0.91$ gives $70.367$ dB, and $x=0.73$--$0.74$ remains near $62$--$63$ dB. Thus the station-design frontier appears locally saturated at this lattice resolution; the next improvement must come from a different station family, a richer validation criterion, or a stronger operator/library model rather than sub-percent phase nudging. We next tested whether more sparse-observed training trajectories improve the shared operator. They do not automatically help. On a fresh test set \{roll $(9,27)$, roll $(27,9)$, roll $(15,29)$, roll $(29,15)$, scale $0.70$, scale $1.50$, scale $0.40$, scale $1.90$\}, the current eight-training-variant row with $8.75\%$ phase $(0.75,0.90)$ reaches $73.381$ dB. Adding four more sparse-observed training variants \{roll $(1,23)$, roll $(23,1)$, scale $0.60$, scale $1.60$\} while keeping the same per-trajectory station budget and solver drops the same fresh-test mean to $70.316$ dB. The support remains exact, but the coefficient vector shifts, especially in the nonlinear and damping terms. Thus the next operator-learning lever is not simply more trajectories; training variants must be selected or weighted so that assimilation bias from scale-extreme trajectories does not distort the shared weak-form coefficients. The isolating controls confirm that the degradation is not caused by one family alone. Adding only the two extra roll variants to the eight-variant training set gives $71.321$ dB on the same fresh test set; adding only the two scale-extreme variants gives $70.131$ dB. Both retain exact support, but both move the coefficients away from the high-PSNR eight-variant estimate. This suggests that the original eight sparse-observed trajectories already form a good coefficient-calibration design for this station phase. Additional trajectories should enter only through validation-selected weights or subset selection, not by unweighted concatenation. We implemented that weighting hook in the runner as \texttt{--train\_variant\_weights}, multiplying each variant's weak-form rows and targets by the square root of its weight before the closed-form solve. Downweighting the four rejected variants improves over unweighted concatenation but still does not beat the eight-variant subset: weights $0.25$, $0.10$, $0.03$, and $0.01$ on the four added variants yield $72.490$, $73.035$, $73.280$, and $73.346$ dB, respectively, on the same fresh test set, versus $73.381$ dB for weight zero. Thus the validation-selected action for these candidates is rejection. The useful research conclusion is that the no-backprop operator learner can support neuromodulatory-style reliability weights, but the first weighted audit says the next gain requires discovering better candidate trajectories or operator features, not softly retaining known harmful variants. We then audited three alternative explanations before changing the operator model. First, a threshold/decoy-pressure sweep around the $8.75\%$ phase $(0.75,0.90)$ frontier used thresholds $10^{-7}$, $2\times10^{-7}$, $5\times10^{-7}$, $10^{-6}$, and $2\times10^{-6}$ with decoy multipliers $5$, $10$, $20$, $50$, and $100$. Low decoy penalties admitted false positives and dropped validation to about $66$ dB, but every exact-support row gave the same $73.568$ dB validation score; applying the validation-selected $10^{-7}$, multiplier-$50$ row to the $12$-variant fresh audit reproduced $73.379$ dB. Thus the frontier is not limited by the sparse threshold once decoys are suppressed. Second, we implemented phase-preserving score-mixed station policies such as \texttt{lattice\_phase75\_90\_gradient01}. These keep the tuned phase lattice as the backbone and replace only $1$--$5\%$ of the station budget with diverse high-gradient, high-variance, or hybrid-score sites. This fairer saliency audit was decisively negative. At the fixed $352$--$358$ station scale, generic score-mixed policies tied to the untuned lattice collapsed to $49.713$ dB or worse, and even the phase-preserving variants degraded monotonically: gradient replacement at $1\%$, $2\%$, $3\%$, and $5\%$ gave $68.468$, $65.134$, $62.046$, and $52.616$ dB; variance replacement gave $55.350$, $51.061$, $46.917$, and $42.942$ dB; hybrid replacement gave $58.324$, $54.872$, $50.958$, and $45.586$ dB. The conclusion is that pointwise saliency is not an adequate station objective for this weak operator learner. The lattice points themselves carry Fourier conditioning, and replacing even a few of them damages the coefficient estimate despite exact support in several rows. The positive improvement came from exact decimation of the phase lattice. Scanning the integer station counts $350$ through $361$ at phase $(0.75,0.90)$ found a new validation winner at $352$ stations, i.e. sensor fraction $352/4096=0.0859375$. This row recovers exact $(7,0,0)$ support and reaches $74.442$ dB on the validation variants, compared with $73.568$ dB for the previous $358$-station row and $73.267$ dB for the complete $19\times19$ grid with $361$ stations. Nearby counts are sharply worse: $350$--$351$ give about $66$ dB, $353$--$354$ give $68.795$--$69.939$ dB, $356$ gives $66.661$ dB, and $359$--$360$ give $70.570$--$72.844$ dB. A local phase refinement at the $352$-station count confirmed $(0.75,0.90)$, with $(0.75,0.91)$ tied by an effectively equivalent rounded station set. On the $12$-variant fresh audit, the validation-selected $352$-station row reaches $73.925$ dB with exact support and a $1.21$ ms sparse solve, improving the previous $358$-station fresh frontier of $73.379$ dB while using fewer observations. This is now the cleanest sparse shallow-water operator-learning result in the manuscript: progress came not from more data, saliency replacement, or threshold tuning, but from validation-selected Fourier-compatible station decimation. We added an exact \texttt{--sensor\_counts} option and widened the decimation sweep to counts $320$--$380$ at the same phase. This exposed an even sharper sparse resonance at $334$ stations, i.e. $334/4096=0.08154296875$ observations per training trajectory. The $334$-station row reaches $75.851$ dB on the validation variants with exact $(7,0,0)$ support, while nearby counts again fluctuate strongly: $328$--$330$ sit near $71$--$72$ dB, $332$ activates false support, $335$ gives $71.424$ dB, $336$ activates two false positives, and the entire $362$--$380$ side-$20$ band stays below $68$ dB except for false-support rows. A fresh $12$-variant audit of the locked $334$-station row reaches $75.841$ dB with exact support and a $1.18$ ms sparse solve. Local phase refinement at count $334$ again selects $(0.75,0.90)$, with $(0.75,0.91)$ tied by the rounded station set. This supersedes the $352$-station checkpoint: the validation-selected operator now improves the broad fresh audit by $2.462$ dB over the previous $358$-station frontier while using $6.7\%$ fewer observations. We also checked whether the same phase contains an even lower-count resonance. A validation sweep over exact station counts $220$--$319$ at phase $(0.75,0.90)$ was negative. The best row in that band is $318$ stations at only $68.349$ dB, and most rows sit near $55$--$66$ dB, with occasional false-support failures such as counts $232$, $235$, $240$, $297$, and $306$. Thus the current sparse optimum is not simply ``as few stations as possible.'' For this Fourier dictionary and weak-form window, the useful resonance appears to start near the high end of the side-$19$ decimation family, with $334$ stations as the current validated minimum-quality sweet spot. The next audit asked whether the $334$-station geometry was limited by the Fourier assimilation basis itself. Holding the training variants, validation variants, station count, phase $(0.75,0.90)$, window, threshold, and decoy pressure fixed, we swept the reconstruction basis from modes $2$ through $6$. Mode $2$ still selected exact support but underfit the observed window and biased the nonlinear coefficient, reaching only $52.464$ dB validation PSNR. The previous mode-$3$ row reached $75.851$ dB validation and $75.841$ dB on the locked $12$-variant fresh audit. Mode $4$ gives a small but clean improvement: it reaches $76.474$ dB on validation, exact $(7,0,0)$ support, and a $1.23$ ms sparse solve; the disjoint $12$-variant audit reaches $76.444$ dB with exact support and a $1.30$ ms solve. Modes $5$ and $6$ regress to $75.486$ and $74.903$ dB, respectively, despite exact support. Repeating the exact-count sweep $320$--$380$ under mode $4$ again selects $334$ stations; count $352$ rises to $74.894$ dB but stays below the $334$-station row, and the side-$20$ band remains weaker or false-support. The current interpretation is therefore a two-axis resonance: the best sparse operator learner is not maximal observation count or maximal basis bandwidth, but the mode-$4$, $334$-station, phase-balanced Fourier geometry. We then revisited the weak temporal projection itself. The committed rows used a window of $7$ frames, stride $2$, and test-mode radius $3$. At the locked mode-$4$, $334$-station geometry, shortening the window is a major coefficient-calibration lever. With stride $2$ and test-mode radius $3$, validation PSNR rises from $76.474$ dB at window $7$ to $76.131$ dB at window $6$, $78.346$ dB at window $5$, $78.244$ dB at window $4$, $79.204$ dB at window $3$, and $79.582$ dB at window $2$, all with exact $(7,0,0)$ support. The projection radius is sharp: at window $3$, radius $2$ collapses to $59.003$ dB and radius $4$ drops to $69.496$ dB; at window $5$, radius $2$ and $4$ give $59.126$ and $70.595$ dB. The lower boundary and stride controls also reject a trivial ``shorter is always better'' rule: window $1$ gives $78.707$ dB, while window $2$ with stride $1$ and $3$ gives $79.189$ and $78.717$ dB. The selected weak setting is therefore window $2$, stride $2$, radius $3$. On the locked $12$-variant fresh audit this reaches $79.509$ dB, exact support, and a $1.26$ ms sparse solve, improving the previous mode-$4$ fresh frontier by $3.065$ dB and the old $358$-station frontier by $6.130$ dB. The inferred vector $(H,\beta,g,f,r,\nu,\mu_h)=(0.99980,0.72348,0.86033,0.58257,0.06370,0.005998,0.004022)$ is now close to the true operator across all seven groups. The lesson is precise: long weak windows were smearing the sparse-assimilated trajectory balance; a short, non-overdense weak window better matches the local truncation and assimilation error scale. Re-sweeping station counts after the weak-window correction shows that this is not a new count-search problem. With mode $4$, window $2$, stride $2$, radius $3$, phase $(0.75,0.90)$, and the same validation variants, counts $300$--$360$ again select $334$ stations at $79.582$ dB. The low-count side remains far below the frontier: the best sub-$320$ count is $318$ at $68.956$ dB. The old $352$-station checkpoint improves from $74.894$ dB to $76.376$ dB under the shorter weak window, and the previous $358$-station phase row rises to $74.607$ dB, but both remain clearly below $334$. Secondary bumps such as count $328$ at $72.560$ dB and count $343$ at $73.083$ dB do not change the ordering. Thus the weak-window correction improves coefficient calibration at fixed geometry, while the station-count resonance itself remains locked. We also ported exact count and phase-station support into the sparse weak-form library-discovery runner, then audited the harder raw and grouped libraries at the locked geometry. This uses the same mode-$4$, $334$-station, phase $(0.75,0.90)$, window-$2$ weak system, but asks the selector to reject nuisance terms rather than assuming the seven physical groups. The raw typed weak solve at threshold $5\times10^{-7}$ and decoy penalty $5$ recovers all $13$ physical columns with zero decoys, active support $(13,0,0)$, in $0.54$ ms. Its grouped typed counterpart recovers all seven shared physical groups with zero decoys, active support $(7,0,0)$, in $0.56$ ms. Raising the threshold to $5\times10^{-6}$ prunes weak true terms, and thresholds $5\times10^{-5}$ or larger over-prune the operator. The forecast values in this single-trajectory discovery runner remain near $20$ dB because the rollout starts from the sparse-assimilated state whose reconstruction PSNR is only about $21$ dB; this row should therefore be read as a support-identifiability result, not as the high-quality multitrajectory forecast frontier above. We then inserted the same raw-vs-grouped choice into the high-quality multitrajectory forecast protocol. This separates support recovery from long-horizon transfer. At the locked mode-$4$, $334$-station, phase $(0.75,0.90)$, window-$2$, stride-$2$, radius-$3$ setting, the raw typed library again recovers exact support $(13,0,0)$ at threshold $5\times10^{-7}$, but its two-variant validation forecast is only $73.676$ dB. The grouped typed operator recovers $(7,0,0)$ and reproduces the $79.582$ dB frontier. Thus the seven-group collapse is not merely a reporting convention. Enforcing the shared physical coefficients $(\beta,g,f,r,\nu)$ across their equation-specific columns is a strong structural regularizer for rollout quality, even when the ungrouped raw support is exactly correct. We also repeated the station and weak-projection controls under the improved window-$2$ setting. A $20$-policy local phase grid with $x\in\{0.65,0.70,0.75,0.80,0.85\}$ and $y\in\{0.80,0.85,0.90,0.95\}$ again selects phase $(0.75,0.90)$ at $79.582$ dB; the nearest strong neighbor is $(0.75,0.80)$ at $79.179$ dB, while many exact-support phases fall into the $60$--$69$ dB range. Phase-preserving score replacement remains decisively negative at the locked count: replacing only $2$--$15\%$ of the $(0.75,0.90)$ lattice by gradient, variance, hybrid, or leverage stations never improves the frontier. The best mixed row is gradient-$2\%$ at $66.530$ dB, and larger gradient replacements introduce decoys or drop below $54$ dB; variance, hybrid, and leverage replacements are similarly weaker. A targeted low-amplitude train-weighting control also loses: weights $(1,1,0.5,1,1,1,2,1)$ on \{base, roll $(7,11)$, scale $1.25$, roll $(13,5)$, roll $(5,17)$, roll $(2,19)$, scale $0.90$, scale $1.10$\} reach $79.153$ dB, below the equal-weight row. Finally, the radius sweep closes the weak-test-mode axis for window $2$: radius $1$ diverges, radius $2$ gives $58.987$ dB, radius $3$ gives $79.582$ dB, radius $4$ gives $69.786$ dB, and radius $5$ gives $65.683$ dB. The selected radius is therefore not arbitrary; it is the unique tested projection scale that balances sparse-assimilation bias and weak-form identifiability. As a no-backprop control-theory follow-up, we added optional sparse-sensor rollout calibration to the multitrajectory runner via \texttt{--theta\_calibration\_steps}. Starting from the weak-form coefficient vector, the routine performs coordinate search over the seven physical parameters and accepts changes that reduce a held-out training-tail loss measured only at the observed sparse station locations. This is biologically and control-theoretically plausible in the sense that it uses forward rollouts and local observation residuals, not reverse-mode differentiation or full-field labels. On the current frontier, however, it is a hard negative result. With steps $1\%$, $0.3\%$, and $0.1\%$, the sparse sensor-tail loss decreases only from $0.2092457$ to $0.2092395$, while validation PSNR collapses from $79.582$ dB to $57.751$ dB. The calibrated vector moves to $(0.99780,0.73364,0.85514,0.58257,0.06281,0.005914,0.003970)$. The conclusion is useful: direct sparse-sensor replay is an overfitting objective for this problem. The weak-form grouped solve generalizes because it optimizes an operator balance, not because it best replays a short sparse observation tail. We then attacked the same bottleneck from the reconstruction and weak-row side. First, temporal smoothing of the assimilated training sequence before weak integration is negative: centered binomial-$3$, binomial-$5$, and box-$3$ smoothing reduce validation to $54.490$, $50.174$, and $51.966$ dB. The smoothed sequences have nearly the same reconstruction PSNR as the raw assimilated fields, but their weak coefficients are biased. Second, concatenating additional weak views is also negative. At the locked mode-$4$, $334$-station, phase $(0.75,0.90)$ geometry, a $50/50$ window-$2$/window-$3$ weak-row mixture reaches only $79.341$ dB, an $80/20$ mixture reaches $79.470$ dB, and a $90.9/9.1$ window-$2$/window-$4$ mixture reaches $79.344$ dB. A lightly weighted radius-$4$ test-mode view also degrades the diagnostic validation run. Thus the selected window-$2$, radius-$3$ weak projection is not merely underdetermined; adding nearby valid projections injects biased rows rather than averaging out the reconstruction error. The next reconstruction audit clarifies the direction. Raising the Fourier assimilation dictionary from mode $4$ to modes $5$, $6$, and $7$ lowers the training-window reconstruction PSNR and drops validation to $76.850$, $75.503$, and $71.546$ dB, even though exact $(7,0,0)$ support is still recovered. In contrast, applying a rectangular low-pass filter to the mode-$4$ assimilated fields before forming the weak rows gives a tiny but clean improvement when the keep radius is $3$. Filter radius $2$ underfits and reaches only $76.787$ dB; radius $4$ is effectively the unfiltered baseline at $79.583$ dB. Radius $3$ reaches $79.587$ dB on validation and $79.519$ dB on the $12$-variant fresh audit with exact support and a $1.19$ ms sparse solve. The coefficient vector is $(0.99980,0.72336,0.86033,0.58258,0.06371,0.005998,0.004022)$. This supersedes the unfiltered same-ridge fresh row at $79.511$ dB and the previous unfiltered frontier at $79.509$ dB, but only by about $0.01$ dB. The scientific value is therefore diagnostic rather than headline: the weak learner is now limited by aliasing and sparse-reconstruction bias at the operator-balance level. We then made this anti-aliasing more operator-specific. The runner now supports \texttt{--weak\_filter\_application} and \texttt{--assim\_spatial\_filter\_shell\_weight}. The selected row keeps the weak target and all linear operator columns on the unfiltered mode-$4$ assimilated sequence, replaces only the true quadratic flux/advection columns by their mode-$3$ low-pass values, and retains the first excluded Fourier shell with weight $0.05$. This is a term-local weak-form filter: it does not smooth the rollout state, does not alter the sparse station observations, and does not use held-out future fields. On the validation variants this \texttt{nonlinear\_terms} row reaches $79.588$ dB with exact $(7,0,0)$ support; on the locked $12$-variant fresh audit it reaches $79.521$ dB, exact support, and a $1.24$ ms sparse solve, with coefficients $(0.99980,0.72336,0.86033,0.58258,0.06371,0.005998,0.004022)$. Filtering nonlinear terms and decoys gives the same validation score; filtering all columns with the same shell reaches only $79.519$ dB fresh. Finally, the same nonlinear-only filter does not make higher assimilation bandwidth safe: mode $5$ and mode $6$ validation runs fall to $76.693$ and $75.301$ dB despite exact support. The interpretation is sharper: the main aliasing source is the quadratic product library, but high-bandwidth assimilated states also bias the linear weak balance and coefficient calibration. The next material step should therefore be an operator-aware station or anti-aliasing objective that keeps the stable mode-$3$ nonlinear products while preserving the useful mode-$4$ linear state information. A decomposition audit then checked whether the filter should act after products are formed or only on one nonlinear physical channel. Product-column filtering is too late: filtering the already-formed nonlinear product columns reaches only about $79.583$ dB validation, essentially the unfiltered row. Momentum-advection-only state filtering is also negative at $79.581$ dB. Mass-flux-only state filtering is the strongest validation row, reaching $79.591$ dB when only the continuity-equation $\beta$ column is formed from the mode-$3$ filtered state. However, this validation gain does not transfer: the hard mass-flux row reaches $79.520$ dB on the $12$-variant fresh audit, and shell weights $0.05$ and $0.10$ also reach only $79.520$ dB. Thus the fresh frontier remains the all-nonlinear state-prefiltered row above. The useful conclusion is methodological: two validation variants can over-rank continuity-specific anti-aliasing, so the next selector must use a broader validation design or a physically derived anti-aliasing criterion rather than a two-trajectory validation score alone. We therefore widened the selector itself before running further filter searches. The broad validation set contains eight additional roll and amplitude variants, \texttt{roll9\_27}, \texttt{roll27\_9}, \texttt{roll15\_29}, \texttt{roll29\_15}, \texttt{scale070}, \texttt{scale150}, \texttt{scale040}, and \texttt{scale190}, while the eight sparse-observed training variants remain fixed. Under this V8 selector, the unfiltered row scores $79.513$ dB, all-column keep-$3$ filtering scores $79.521$ dB, all nonlinear-state filtering scores $79.523$ dB, the all-nonlinear shell-$0.05$ row scores $79.523$ dB, mass-flux-only hard filtering scores $79.522$ dB, and mass-flux-only shell-$0.05$ filtering scores $79.522$ dB, all with exact $(7,0,0)$ support. The broader selector therefore chooses the same all-nonlinear shell-$0.05$ rule that transferred best to the $12$-variant fresh audit, rather than the continuity-only row that won the narrow two-variant validation. This is the current robust selection rule: keep the mode-$4$ assimilated state for targets and linear operators, form all true nonlinear state products from the mode-$3$ state with a $0.05$ first-shell taper, and validate across both phase rolls and amplitude extremes. Two follow-up geometry audits closed the obvious remaining local axes. Re-sweeping exact station counts $318,328,334,343,352,358$ under the V8 selector and the nonlinear shell rule again selects $334$ stations. The tested rows all recover exact support, but coefficient calibration is sharply count-dependent: the V8 means are $68.898$, $71.748$, $79.523$, $73.055$, $75.713$, and $74.287$ dB. A nearby phase grid at count $334$ is even sharper. The nine policies \texttt{lattice\_phase70\_85} through \texttt{lattice\_phase80\_95} all recover exact support, but only \texttt{lattice\_phase75\_90} reaches the frontier. The other phase rows range from $58.980$ to $68.973$ dB. Thus the station objective is not ``recover the seven groups''; it is to preserve the Fourier-compatible weak-row geometry that calibrates the nonlinear and damping coefficients. We also tested whether the weaker high-amplitude validation rows could be fixed by adding amplitude-extreme trajectories to the training weak solve. On a disjoint V8b holdout, the fixed eight-trajectory training set scores $79.536$ dB. Adding \texttt{scale070} and \texttt{scale150} with equal weights drops to $76.522$ dB; giving those added variants only $0.1$ weight still drops to $79.317$ dB. A new whole-trajectory \texttt{--train\_block\_normalization} control was added to test block-level scaling without row-wise physics destruction. Target-RMS block normalization on the same augmented set activates all $11$ decoys and drops to $65.019$ dB, while the dense true-support coefficient row is still only $77.358$ dB. The result is negative but useful: broad amplitude coverage is not automatically a regularizer for this sparse weak system. The next structural lever should be an operator-domain anti-aliasing or station-design criterion, not more amplitude variants or block/row rescaling. The first useful station-design criterion came from the operator domain rather than from reconstruction or row conditioning. We added \texttt{--operator\_reference\_mode dense\_train}, which fits a dense weak-form reference operator on the training window only and scores each sparse station geometry by the relative distance between its sparse-derived $\theta$ and this dense-training $\theta$. The dense-training reference for the current eight-trajectory set is $(1.00037,0.71970,0.86033,0.58005,0.06454,0.006010,0.004011)$, within $4.59\times10^{-4}$ relative error of the known simulator coefficients. This training-only objective selects \texttt{lattice\_phase75\_90} in the nine-phase grid and selects the same count-$334$, phase-$75\_90$ row in a local $3\times3$ count/phase grid. The selected row has operator-reference error $0.00283$ and V8 PSNR $79.523$ dB; the next closest local candidate is count $352$, phase $75\_90$ with error $0.01088$ and $75.713$ dB. This is not a new quality frontier, but it is a more principled path to station placement: design sparse observation geometry to reproduce the dense training operator, then evaluate transfer on disjoint futures. We then stress-tested that rule as an actual selector over a broader $30$-candidate grid: station counts $318,328,334,343,352,358$ crossed with phases \texttt{70\_90}, \texttt{75\_85}, \texttt{75\_90}, \texttt{75\_95}, and \texttt{80\_90}. On the V8 holdout, both held-out PSNR and the training-only dense-operator objective select count $334$, phase \texttt{75\_90}, with $79.523$ dB and operator-reference error $0.00283$. The top held-out alternatives are count $352$, phase \texttt{75\_90} at $75.713$ dB and count $358$, phase \texttt{75\_90} at $74.287$ dB. Repeating the same grid on a disjoint V8b holdout gives the same top row: $79.536$ dB for count $334$, phase \texttt{75\_90}; the next held-out rows are count $352$, phase \texttt{75\_90} at $75.834$ dB and count $358$, phase \texttt{75\_90} at $74.347$ dB. The caveat is that operator-reference error is a successful top-$1$ selector here, not yet a calibrated total ordering: for example, count $358$, phase \texttt{75\_85} has the second-lowest reference error but only about $68$ dB. The next station-design objective should therefore keep dense-operator matching as the primary constraint, but add a training-only stability term that penalizes geometries whose coefficient match is fragile under small phase, count, or trajectory perturbations. A follow-up diagnostic showed that the failure mode is not local instability but coefficient scaling. The original dense-operator score is an absolute relative L2 distance, so it is dominated by the large coefficients $(H,\beta,g,f)$ and can underweight small but rollout-sensitive coefficients such as damping and viscosity. We therefore added \texttt{--operator\_reference\_metric} with \texttt{absolute\_l2}, \texttt{relative\_l2}, and \texttt{relative\_linf} options. The \texttt{relative\_linf} metric scores the maximum coefficient-wise relative error to the dense training operator. On the same V8 $30$-candidate grid, \texttt{relative\_linf} again selects count $334$, phase \texttt{75\_90}, with score $0.01294$ and $79.523$ dB. The row count $358$, phase \texttt{75\_85}, which was second-best under absolute L2, is demoted to score $0.04446$ because its small-coefficient errors are large. This makes the selector more physically balanced: it still does not perfectly rank all held-out PSNR values, but it removes the most obvious scale artifact in dense-operator matching. Finally, a broader relative-\(\ell_\infty\) selector search crossed eight counts, $300,318,328,334,343,352,358,372$, with the full $3\times3$ phase neighborhood from \texttt{70\_85} through \texttt{80\_95}. Across all $72$ rows, both held-out PSNR and the balanced dense-operator score again select count $334$, phase \texttt{75\_90}, with $79.523$ dB and score $0.01294$. The best non-frontier held-out row is count $352$, phase \texttt{75\_90} at $75.713$ dB; count $372$ never exceeds $67.353$ dB. This closes the obvious count/phase lattice search. The next lever should be the weak-form estimator itself: multi-window test functions, coefficient-balanced row construction, or model-bias correction for the nonlinear product columns. We then held the selected station geometry fixed and tested whether concatenated weak views could reduce estimator bias. They did not. The baseline single view, window $2$ and test radius $3$, remains $79.523$ dB with balanced operator score $0.01294$. Equal-weight windows $(1,2,3)$ drop to $79.334$ dB; equal-weight windows $(2,3,4)$ drop to $78.918$ dB; a conservative $90/10$ mix of windows $(2,3)$ still drops to $79.476$ dB. Spatial test-mode concatenation is worse: modes $(2,3)$ drop to $71.673$ dB and modes $(3,4)$ drop to $74.947$ dB. This also validates the balanced score: the $(1,2,3)$ row has a smaller absolute operator distance than the baseline, but worse balanced relative error and worse rollout. Thus the useful weak-form view is now sharply identified as window $2$, test radius $3$; the remaining bias is not fixed by naive multi-scale row concatenation. Finally, we tested whether the remaining gap is fundamental or merely a calibration-objective problem. A diagnostic forward-only coordinate search was run from the current sparse coefficient vector, using full-field training-tail rollout residuals over the eight training variants and no reverse-mode differentiation. This is not a sparse-only claim, because the calibration objective uses full-field training tails; it is a boundary diagnostic. The result is decisive: the sparse vector $(0.99980,0.72336,0.86033,0.58258,0.06371,0.005998,0.004022)$ gives $79.523$ dB on the disjoint V8 holdout, while the tracked full-field calibration runner \texttt{nonlinear\_shallow\_water\_theta\_calibration\_diagnostic.py} reaches $(1.00000,0.72076,0.85999,0.58002,0.06489,0.006001,0.004000)$, $104.088$ dB on V8, and $104.319$ dB on the disjoint V8b holdout. The accepted moves mainly correct $\beta$, $f$, $r$, $H$, and $\mu_h$ toward the simulator values. In contrast, rolling out the dense-training weak-form reference vector, although closer to the simulator coefficients, gives only $78.583$ dB. Therefore the $79.5$ dB frontier is not a support, station-count, or weak-window ceiling; it is a coefficient-calibration objective ceiling. The follow-up calibration audit isolates what information is missing. If the calibration starts from the exact full training-tail state but observes only the $334$ sparse station values over the tail, the same forward-only coordinate search reaches the same vector and the same V8/V8b values, $104.088/104.319$ dB. Sparse station values therefore contain enough coefficient information once the hidden state at the calibration start is correct. The failure is the state anchor: using the sparse-assimilated tail-start state makes even a full-field tail loss collapse to $55.875$ dB, and optimizing either the sparse station tail, the assimilated pseudo-field tail, or the weak residual of the assimilated tail drives the coefficients in the wrong direction. Dense-training weak residual and dense-operator relative-$\ell_\infty$ objectives are also not sufficient by themselves: they reduce their training losses but reach only $78.997$ and $79.294$ dB on V8. We then tested sparse or partial state-anchor repairs. A low-mode forward-sensitivity correction fitted from station history is too weak: mode-$2$, mode-$3$, and mode-$4$ anchors improve full-tail PSNR by only $0.00029$, $0.00043$, and $0.00084$ dB, respectively, and the calibrated mode-$2$ row still collapses to $63.281$ dB. A stronger but more privileged dynamics anchor, obtained by propagating the known full initial training state to the calibration split with the sparse OSNR operator, gives a much better training-tail state ($75.778$ dB tail PSNR and station loss $5.98\times10^{-7}$). Unconstrained station-tail calibration from this anchor still overfits, dropping to $73.021$ dB after $39$ accepted moves, but a one-accepted-move trust-region update raises only the damping coefficient, $r:0.0637088\mapsto0.0643459$, and improves V8/V8b to $80.123/80.139$ dB. A two-move variant immediately accepts an overlarge Coriolis correction and falls to $76.831$ dB. Thus the first legitimate direction beyond the $79.5$ dB row is not blind coordinate search; it is trust-region, state-aware station calibration. The next selector audit made that trust-region rule explicit in the tracked diagnostic runner. The new one-step selector enumerates all single-coordinate candidates at the configured step sizes, ranks them using training-side objectives only, and evaluates V8/V8b only after the selected move is fixed. Station-tail loss alone is a negative control: it selects a depth decrease, $H:0.9998017\mapsto0.9988019$, because that gives the largest station-tail improvement, but the held-out scores collapse to $74.987/74.999$ dB. Adding a dense training-tail weak-residual consistency gate changes the selected move. Among candidates that improve both the model-anchor station tail and the dense training weak residual, the selector chooses the small Coriolis correction $f:0.5825801\mapsto0.5808324$; only after that training-only selection do the disjoint diagnostics evaluate to \textbf{$82.194/82.221$ dB} on V8/V8b. This is a larger lift than the previous damping-only trust move, but it is still diagnostic rather than sparse-only, because the auxiliary gate uses dense training-tail weak rows. Multi-accept variants do not improve the result: accepting a subsequent viscosity move gives $82.050/82.079$ dB, and an additional dense-operator relative-$\ell_2$ gate still keeps the best state at the first accepted $f$ move. The useful conclusion is sharper: sparse station replay supplies candidate moves but cannot select them safely by itself; a training-side operator-consistency gate can reject destructive station overfits and select a real coefficient correction. The next publishable step is to replace the dense weak-residual gate with a station-observable or assimilated operator-consistency surrogate while preserving the one-move trust-region discipline. That replacement attempt is now also informative. We added station-observable selector objectives based on sparse-assimilated weak residuals, station finite-difference RHS residuals, one-step station-increment replay, held-out station splits, and separate higher-mode gate reconstructions. None is a safe substitute for the dense weak gate. The sparse weak gate selects the same destructive $H$ move and gives $74.987/74.999$ dB; station-RHS and station-one-step primaries select $\beta:0.7233627\mapsto0.7161291$ and slightly reduce V8/V8b to $79.489/79.509$ dB; held-out station one-step replay selects $g:0.8603310\mapsto0.8517277$ and collapses to $55.525/55.525$ dB; higher-mode sparse weak gates at modes $5$ and $6$, with and without temporal smoothing, still admit the destructive $H$ move. A clean sparse weak-solve hyperparameter audit over shell weights $0,0.025,0.075,0.10$ and sensor ridges $10^{-5},10^{-4}$ also fails to move the frontier, staying at $79.520$--$79.523$ dB. The only positive replacement so far is structural rather than learned: restrict the selector to momentum coefficients $(f,r,\nu)$ and to small trust steps $(0.003,0.001)$. With the same model-anchor station primary, this selects the same $f:0.5825801\mapsto0.5808324$ move and reaches $82.194/82.221$ dB without dense weak rows; on two fresh eight-variant audit lists it moves $79.521/79.524$ dB to $82.191/82.194$ dB. This is a cleaner diagnostic than the dense-gated selector, but it is still not a sparse-only claim because the selector primary uses the model-propagated full initial training state. When the primary is changed to the fully sparse assimilated station tail, the same small-trust momentum selector chooses $\nu:0.0059976\mapsto0.0059796$ and drops to $78.351/78.360$ dB. The next real problem is therefore state anchoring: station observations contain the coefficient signal, but the current sparse-assimilated calibration-start state distorts the selector enough that station-local objectives prefer wrong coefficient directions. The state-anchor follow-up gives the first clean sparse calibration lift. Instead of using the privileged full initial training state, we propagate only the sparse-assimilated initial frame to the calibration split with the learned sparse OSNR operator, then run the same one-move small-trust momentum selector on station-tail loss. This \emph{sparse-model anchor} is still a low-PSNR field in full space (mean start/tail PSNR $21.069/20.386$ dB over the training variants), but it is dynamically consistent with the learned operator and improves the station-tail selector geometry. With primary objective \texttt{sparse\_model\_station\_tail}, coordinate subset $(f,r,\nu)$, and steps $(0.003,0.001)$, the selector chooses $f:0.5825801\mapsto0.5808324$ and moves V8/V8b from $79.523/79.536$ dB to \textbf{$82.194/82.221$ dB}, without dense weak rows, without full-field calibration tails, and without the full-initial-state model anchor. The same selected move transfers on two fresh eight-variant audit lists, $79.521/79.524\mapsto82.191/82.194$ dB. The constraints are sharp: allowing the large $0.01$ step without validation oversteps to $f=0.5767543$ and falls to $76.540/76.547$ dB; allowing a second small momentum move falls to $80.905/80.925$ dB; allowing all coordinates even at small trust selects $\beta:0.7233627\mapsto0.7255328$ and falls to $79.213/79.248$ dB. We then added a training-variant split to the sparse-model station objective. The selector can now use \texttt{sparse\_model\_station\_tail\_fit} as the primary loss and require improvement on \texttt{sparse\_model\_station\_tail\_val}. This removes the hand-coded step scale within the momentum subspace: with candidate steps $(0.01,0.003,0.001)$, the validation split rejects the destructive $0.01$ overstep and selects the same small $f:0.5825801\mapsto0.5808324$ move, preserving $82.194/82.221$ dB. A second split-validated refinement from that point, using only $f$ and smaller steps $(0.001,0.0003,0.0001,0.00003)$, accepts $f:0.5808324\mapsto0.5806582$ and improves V8/V8b to \textbf{$82.278/82.305$ dB}; a fresh disjoint two-list audit gives $82.274/82.278$ dB. A subsequent $(f,r,\nu)$ pass finds no eligible move. However, the split does \emph{not} learn the coordinate mask: all-coordinate split validation at small trust selects $H$ and collapses to $66.444/66.443$ dB, all-coordinate large-trust split validation selects $\beta$ and gives $78.174/78.270$ dB, and adding a sparse weak-residual gate selects $g$ and collapses to $66.058/66.061$ dB. The publishable statement is therefore narrow but stronger than before: a sparse station-derived dynamic anchor, a physics-motivated momentum coordinate mask, and a training-split trust rule give a reproducible $+2.75$ dB held-out lift over the $334$-station exact-support frontier. The next missing piece is a station-observable coordinate-confidence rule, likely based on operator-block sensitivity or adjoint/Fisher geometry, that rejects the compensatory $H,\beta,g$ moves without using dense labels or held-out futures. \subsection{Sparse governing-equation discovery for nonlinear shallow water} The preceding experiment assumes that the correct seven-column operator family is known. The more ambitious physics-learning problem is to discover the governing equation itself from a larger nonlinear library. We implemented \texttt{apps\_industrial\_breakthrough/nonlinear\_shallow\_water\_library\_discovery.py}, which expands the residual library to $24$ candidate columns: the $13$ true equation-specific terms for height and both velocity channels, plus $11$ decoys including raw fields, quadratic field products, and misplaced height-gradient terms. OSNR applies a sequential thresholded least-squares solve on only $0.25\%$ of the training residual rows. The discovered support is then collapsed back into the shared physical vector $(H,\beta,g,f,r,\nu,\mu_h)$ and used for held-out future forecasting. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Setting & Support $(TP,FP,FN)$ & Future PSNR & Discovery time & Adam support & Adam PSNR & Adam time \\ \midrule Clean, STLS threshold $5\cdot10^{-4}$ & $(13,0,0)$ & $72.7835$ dB & $1.0432$ ms & $(10,11,3)$ & $43.8385$ dB & $1587.90$ ms \\ $0.1\%$ noise, band-$18$ & $(13,0,0)$ & $69.2189$ dB & $1.4878$ ms & $(10,10,3)$ & $33.8352$ dB & $1632.81$ ms \\ $0.2\%$ noise, band-$16$ & $(12,0,1)$ & $64.5354$ dB & $1.6853$ ms & $(12,0,1)$ & $66.7134$ dB & $1634.47$ ms \\ Oracle true-library LS, clean & n/a & $72.6779$ dB & $0.2817$ ms & n/a & n/a & n/a \\ \bottomrule \end{tabular} \caption{Sparse nonlinear shallow-water governing-equation discovery from a $24$-term overcomplete library. OSNR recovers the exact clean PDE support and remains robust at $0.1\%$ observation noise. The Adam library baseline uses the same anchor rows and $20{,}000$ optimization steps with an $\ell_1$ penalty.} \label{tab:nonlinear-shallow-water-library-discovery} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/nonlinear_shallow_water_library_discovery_outputs_thr0005_adam20k/nonlinear_shallow_water_library_discovery.png} \caption{Sparse governing-equation discovery for the coupled nonlinear shallow-water core. The discovered PDE recovers the clean future at $72.7835$ dB after selecting all true terms and no decoys from the overcomplete library.} \label{fig:nonlinear-shallow-water-library-discovery} \end{figure} This is the first result in the project that begins to look like a genuine PINN/SINDy-class breakthrough rather than only a fast solver. In the clean case, OSNR selects every true nonlinear PDE term and rejects every decoy, then slightly outperforms the oracle true-library forecast. The comparable Adam library run is not just slower; after $20{,}000$ gradient steps it keeps $11$ false-positive decoys, misses $3$ true terms, and loses almost $29$ dB of forecast quality. The wall-clock ratio for discovery is about $1522\times$ in favor of OSNR. At $0.1\%$ observation noise, the same sparse support is still recovered exactly and the speed ratio remains above $1000\times$. At $0.2\%$ noise, spectral denoising preserves zero false positives but one weak damping term drops below threshold; the optimized Adam library catches up in forecast quality only after paying the full $1.6$ s optimization cost. This identifies the next hard technical layer: noise-aware thresholding or group sparsity for weak physical terms, not larger neural networks. \subsection{Canonical PDE discovery: Burgers and Kuramoto--Sivashinsky} To reduce the risk that the shallow-water result is viewed as a repository-specific construction, we added \texttt{apps\_industrial\_breakthrough/canonical\_pde\_discovery\_benchmark.py}. It evaluates two standard equation-discovery controls: viscous Burgers, \[ u_t=-u u_x+\nu u_{xx}, \] and the chaotic Kuramoto--Sivashinsky equation, \[ u_t=-u u_x-u_{xx}-u_{xxxx}. \] Both are discovered from the same $10$-term library \[ \{u,u^2,u_x,u u_x,u^2u_x,u_{xx},u u_{xx},u^3,u_{xxx},u_{xxxx}\}, \] using only $128$ sampled residual rows. The Kuramoto--Sivashinsky trajectory is generated with the standard ETDRK4 spectral integrator; the discovery stage is independent of that generator and sees only the sampled field values. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Equation & OSNR support & OSNR PSNR & OSNR time & Adam support & Adam time \\ \midrule Burgers, clean & $(2,0,0)$ & $43.1187$ dB & $0.3690$ ms & $(2,0,0)$ & $1265.03$ ms \\ Kuramoto--Sivashinsky, clean & $(3,0,0)$ & $12.2751$ dB & $0.0972$ ms & $(3,0,0)$ & $1244.63$ ms \\ \bottomrule \end{tabular} \caption{Canonical PDE discovery controls. Support is reported as $(TP,FP,FN)$ against the known governing equation. Adam uses the same sampled residual rows and $20{,}000$ $\ell_1$-regularized optimization steps.} \label{tab:canonical-pde-discovery} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/canonical_pde_discovery_outputs_clean_thr01/canonical_pde_discovery.png} \caption{Canonical Burgers and Kuramoto--Sivashinsky discovery from a shared overcomplete library. Both OSNR and Adam recover the clean support, but OSNR does it via a millisecond-scale sparse solve rather than a long gradient-optimization loop.} \label{fig:canonical-pde-discovery} \end{figure} The canonical control confirms that the sparse operator-discovery mechanism is not confined to the shallow-water generator. On Burgers, OSNR recovers exactly $\{u u_x,u_{xx}\}$ and is about $3428\times$ faster than the Adam library optimizer. On Kuramoto--Sivashinsky, OSNR recovers exactly $\{u u_x,u_{xx},u_{xxxx}\}$ and is about $12809\times$ faster. The chaotic KS forecast PSNR is naturally low over the held-out horizon because small coefficient and phase errors amplify quickly; for this control, support recovery and coefficient recovery are the meaningful scientific-discovery metrics. At $0.1\%$ direct observation noise, both OSNR and Adam pick decoys under simple pointwise derivative regression, which confirms that the next publishable robustness layer must be weak-form or group-sparse denoised discovery rather than more gradient steps. \subsection{Weak-form canonical PDE discovery under observation noise} The pointwise canonical experiment exposes the correct failure mode: differentiating noisy data directly creates spurious high-frequency library columns. We therefore implemented the weak-form variant in \texttt{apps\_industrial\_breakthrough/canonical\_pde\_weakform\_discovery.py}. Instead of regressing $u_t$ at individual grid points, OSNR integrates the PDE over temporal windows and projects the resulting balance onto low-frequency spatial Fourier test functions: \[ u(t_b)-u(t_a)=\int_{t_a}^{t_b}\Theta(u(t))\,\xi\,dt. \] This is the operator-spline analogue of weak-form PDE discovery: the test functions absorb observation noise before sparse regression sees the library. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Equation/noise & OSNR support & OSNR PSNR & OSNR time & Adam support & Adam PSNR & Adam time \\ \midrule Burgers, $0.1\%$ & $(2,0,0)$ & $49.8283$ dB & $0.4965$ ms & $(2,0,0)$ & $47.9135$ dB & $1560.30$ ms \\ Burgers, $0.5\%$ & $(2,0,0)$ & $49.3129$ dB & $0.2580$ ms & $(2,0,0)$ & $40.7179$ dB & $1537.55$ ms \\ Burgers, $1.0\%$ & $(2,0,0)$ & $48.8199$ dB & $0.2103$ ms & $(2,0,0)$ & $45.2210$ dB & $1539.72$ ms \\ KS, $0.1\%$ & $(3,0,0)$ & $12.6155$ dB & $0.2870$ ms & $(3,0,0)$ & $12.7014$ dB & $1549.60$ ms \\ KS, $0.5\%$ & $(3,0,0)$ & $12.5214$ dB & $0.3026$ ms & $(3,0,0)$ & $11.9874$ dB & $1560.78$ ms \\ KS, $1.0\%$ & $(3,0,0)$ & $12.6433$ dB & $0.3148$ ms & $(3,0,0)$ & $11.7123$ dB & $1529.31$ ms \\ \bottomrule \end{tabular} \caption{Weak-form canonical PDE discovery under observation noise. The support tuple is $(TP,FP,FN)$. Both methods use the same weak rows and library, while Adam uses $20{,}000$ $\ell_1$-regularized optimization steps.} \label{tab:weak-canonical-pde-discovery} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/canonical_pde_weakform_outputs_thr002/canonical_pde_weakform_discovery.png} \caption{Weak-form Burgers and Kuramoto--Sivashinsky discovery under noisy observations. The weak operator rows recover the correct governing support through $1\%$ noise while avoiding the pointwise derivative decoys observed in the direct regression control.} \label{fig:weak-canonical-pde-discovery} \end{figure} This is the strongest canonical scientific-ML result so far. The weak-form OSNR solver recovers the exact Burgers and KS support at every tested noise level up to $1\%$. On Burgers, it is also materially more accurate than Adam in forecast quality, improving the $0.5\%$ noise row by $8.5950$ dB. On KS, both methods recover the same support, but OSNR reaches the solution roughly $4{,}858\times$ to $5{,}400\times$ faster for the threshold-$0.002$ profile. This converts the earlier ``mostly speed'' canonical result into a robustness result: integral operator rows eliminate noisy derivative decoys while preserving millisecond-scale discovery. \subsection{Weak-form two-dimensional Navier--Stokes vorticity discovery} The next CFD-facing step is a genuinely two-dimensional incompressible flow operator rather than a scalar one-dimensional PDE. We implemented \texttt{apps\_industrial\_breakthrough/navier\_stokes\_weakform\_discovery.py}, which generates periodic vorticity trajectories and discovers the vorticity equation \[ \omega_t = -u\omega_x-v\omega_y+\nu\Delta\omega, \qquad u=\psi_y,\quad v=-\psi_x,\quad -\Delta\psi=\omega. \] The discovery stage is not told the two-term equation. It sees a $12$-term library containing the advective term, the Laplacian, and ten decoys built from raw vorticity, velocity, first derivatives, and nonlinear products. As in the canonical weak-form experiment, OSNR integrates over time windows and projects the balance onto low-frequency two-dimensional Fourier test functions, \[ \omega(t_b)-\omega(t_a)=\int_{t_a}^{t_b}\Theta(\omega(t),u(t),v(t))\,dt, \] then applies a scaled sequential thresholded solve. The Adam control optimizes the same weak rows for $20{,}000$ $\ell_1$-regularized steps. A first high-resolution attempt at $160^2$ with the coarse timestep became numerically unstable, so the retained scaled run tightens the timestep and reference substepping rather than hiding the CFL boundary. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Grid/noise & OSNR support & OSNR coefficients $(c_{\rm adv},\nu)$ & OSNR PSNR & OSNR time & Adam support & Adam time \\ \midrule $96^2$, clean & $(2,0,0)$ & $(0.9999865,\;0.0015000)$ & $135.3978$ dB & $0.6477$ ms & $(2,11,0)$ & $1958.15$ ms \\ $96^2$, $0.1\%$ & $(2,0,0)$ & $(0.9999592,\;0.0014999)$ & $121.0461$ dB & $0.5487$ ms & $(2,11,0)$ & $2013.22$ ms \\ $96^2$, $0.5\%$ & $(2,0,0)$ & $(1.0007806,\;0.0015010)$ & $99.6485$ dB & $0.3822$ ms & $(2,8,0)$ & $1966.74$ ms \\ $128^2$, clean & $(2,0,0)$ & $(0.9999944,\;0.0015000)$ & $146.9002$ dB & $0.7231$ ms & $(2,11,0)$ & $2322.75$ ms \\ $128^2$, $0.1\%$ & $(2,0,0)$ & $(1.0000714,\;0.0015000)$ & $126.0550$ dB & $0.4781$ ms & $(2,11,0)$ & $2324.12$ ms \\ \bottomrule \end{tabular} \caption{Weak-form two-dimensional Navier--Stokes vorticity discovery. The support tuple is $(TP,FP,FN)$ relative to the two true terms $\{-u\omega_x-v\omega_y,\Delta\omega\}$. Adam is the same weak-library regression optimized by backpropagation, not a full neural Navier--Stokes model.} \label{tab:navier-stokes-weakform-discovery} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/navier_stokes_weakform_outputs_128_stable/navier_stokes_weakform_discovery.png} \caption{Scaled $128^2$ Navier--Stokes weak-form discovery. OSNR identifies the exact advection--diffusion vorticity operator and forecasts the held-out future from the recovered coefficients, while the gradient-optimized sparse regression admits many decoys.} \label{fig:navier-stokes-weakform-discovery} \end{figure} This result is the first high-impact two-dimensional CFD discovery benchmark in the repository. It is still a controlled periodic vorticity system, not a direct DeepMind weather-model comparison. The important claim is narrower and stronger: for a known candidate library and noisy observations, weak-form OSNR recovers the exact incompressible Navier--Stokes vorticity support and coefficients at $96^2$ and $128^2$ resolution, remains stable through $0.5\%$ noise in the $96^2$ run, and solves the sparse operator identification in less than a millisecond. The Adam control uses the same rows and library but remains thousands of times slower and selects many decoy terms. The next step is therefore to move from periodic vorticity discovery to partial-observation assimilation and forced/stochastic Navier--Stokes, where the sparse innovation machinery can be tested on genuinely unknown forcing rather than only coefficient recovery. \subsection{Forced Navier--Stokes sparse innovation assimilation} We then tested the more realistic assimilation problem in \texttt{apps\_industrial\_breakthrough/navier\_stokes\_sparse\_forcing\_assimilation.py}. The governing operator is assumed known, but the trajectory is driven by hidden sparse spatiotemporal forcing: \[ \omega_t = -u\omega_x-v\omega_y+\nu\Delta\omega+f(t,x,y), \qquad f(t,x,y)=\sum_{r=1}^R a_r\,\varphi_t(t-\tau_r)\varphi_x(x-x_r,y-y_r). \] This is closer to weather and flow data assimilation than coefficient discovery: the unknowns are localized forcing events, not just scalar PDE coefficients. OSNR first denoises the observed trajectory spectrally, applies the known Navier--Stokes operator to form the innovation residual \[ \widehat f(t+\tfrac12)=\frac{\omega(t+\Delta t)-\omega(t)}{\Delta t} -\left[-u\omega_x-v\omega_y+\nu\Delta\omega\right]_{t+1/2}, \] then runs a weak three-dimensional matched atom detector: the residual is convolved with the separable spatial--temporal Gaussian test function associated with the forcing atom before non-maximum suppression. A final small least-squares amplitude debiasing step fits the active atoms to the residual. The comparison baselines are an unforced rollout and a smooth low-pass residual forcing field. \begin{table}[h] \centering \small \begin{tabular}{lcccccc} \toprule Run & Events & Noise & Event error & Recall@3 & Traj. PSNR & Assimilation time \\ \midrule $96^2\times121$ & $48$ & $0.2\%$ & $0.7416$ & $97.92\%$ & $81.9812$ dB & $225.76$ ms \\ $128^2\times161$ & $96$ & $0.2\%$ & $1.2556$ & $94.79\%$ & $83.0769$ dB & $851.00$ ms \\ $128^2\times161$ & $96$ & $0.5\%$ & $2.6500$ & $87.50\%$ & $79.5815$ dB & $883.36$ ms \\ $128^2\times161$ & $96$ & $1.0\%$ & $5.2118$ & $73.96\%$ & $74.5790$ dB & $905.40$ ms \\ $192^2\times201$ & $160$ & $1.0\%$ & $5.0134$ & $78.75\%$ & $79.1944$ dB & $4086.68$ ms \\ $192^2\times201$ & $240$ & $1.0\%$ & $3.8576$ & $82.50\%$ & $77.9266$ dB & $6200.45$ ms \\ $192^2\times201$ & $240$ & $2.0\%$ & $4.6710$ & $77.50\%$ & $70.9027$ dB & $6158.36$ ms \\ $192^2\times201$ & $240$ & $3.0\%$ & $4.8846$ & $77.08\%$ & $67.3531$ dB & $6367.19$ ms \\ $192^2\times201$ & $240$ & $5.0\%$ & $8.1801$ & $60.83\%$ & $61.7619$ dB & $6334.94$ ms \\ \bottomrule \end{tabular} \caption{Forced two-dimensional Navier--Stokes sparse innovation assimilation. Event error is measured in joint $(t,x,y)$ grid units against the injected forcing centers. The weak three-dimensional atom detector keeps the sparse residual useful through $5.0\%$ observation noise and improves trajectory reconstruction over both unforced and low-pass residual rollouts.} \label{tab:navier-stokes-sparse-forcing} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/navier_stokes_sparse_forcing_outputs_192_events240_noise030_smooth/navier_stokes_sparse_forcing.png} \caption{Forced Navier--Stokes sparse innovation assimilation on the $192^2\times201$ stress run with $240$ hidden events at $3.0\%$ observation noise. The recovered weak-form sparse atoms preserve the assimilated trajectory substantially better than an unforced model and better than a smooth low-pass residual forcing field.} \label{fig:navier-stokes-sparse-forcing} \end{figure} The $128^2$ dense-event case recovers $94.79\%$ of the hidden forcing events within three grid units at $0.2\%$ noise and reaches $83.0769$ dB trajectory PSNR, compared with $77.8847$ dB for low-pass forcing and $71.1900$ dB for the unforced operator. After replacing the point detector with the weak three-dimensional matched atom score, the same dense case remains useful at $0.5\%$ and $1.0\%$ noise: at $1.0\%$ noise it recovers $73.96\%$ of events within three grid units and reaches $74.5790$ dB, compared with $69.8860$ dB for the low-pass residual and $67.8846$ dB for the unforced operator. The larger $192^2\times201$ stress run with $240$ hidden events at $1.0\%$ noise recovers $82.50\%$ of events within three grid units and reaches $77.9266$ dB, compared with $72.5998$ dB for the low-pass residual and $70.5923$ dB for the unforced model. With scale-adjusted weak atom smoothing, the same $240$-event stress run remains ahead of the low-pass residual at $2.0\%$, $3.0\%$, and $5.0\%$ observation noise. At $3.0\%$ noise, it recovers $77.08\%$ of forcing events and improves trajectory quality by $3.8774$ dB over low-pass; at $5.0\%$ noise, it still recovers $60.83\%$ of events and keeps a $2.0572$ dB trajectory advantage. This is a meaningful step beyond coefficient discovery: sparse OSNR innovations can assimilate unknown localized forcing in a nonlinear two-dimensional flow under noisy observations. The remaining bottleneck is now external benchmark standardization and heavy-overlap amplitude calibration, not basic sparse forcing recovery. To compare against a trained coordinate-field alternative, we added \texttt{apps\_industrial\_breakthrough/navier\_stokes\_neural\_forcing\_baseline.py}. The neural baseline receives the same innovation residual as OSNR and fits a Fourier-feature MLP $g_\theta(t,x,y)$ with AdamW, using a sample distribution biased toward high residual magnitude so that sparse events are not hidden by uniform sampling. The learned forcing is then rolled through the same Navier--Stokes solver. On the $128^2\times161$ case with $96$ hidden events and $3.0\%$ observation noise, a $5$-layer, $128$-hidden-unit Fourier MLP trained for $5{,}000$ steps reaches only $42.9701$ dB trajectory PSNR. OSNR reaches $64.6633$ dB from the same residual, while the low-pass residual baseline reaches $60.8749$ dB. The recovery step takes $0.9165$ s for OSNR versus $52.1354$ s for the neural training loop, a measured $56.9\times$ speed advantage before rollout. The conclusion is not that this small MLP is a definitive neural SOTA baseline; rather, it isolates the key mechanism: dense coordinate-field training smooths or misallocates sparse innovations, while the operator-sparse residual directly preserves the hidden forcing events. \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/navier_stokes_neural_forcing_outputs_128_noise030_steps5000/navier_stokes_neural_forcing.png} \caption{Forced Navier--Stokes sparse forcing recovery against a trained Fourier-feature neural residual field. The neural field is trained directly on the same residual observations but remains much less accurate in the downstream flow rollout.} \label{fig:navier-stokes-neural-forcing-baseline} \end{figure} Finally, we converted the high-noise forced-flow result into a replicated stress suite in \texttt{apps\_industrial\_breakthrough/navier\_stokes\_high\_noise\_suite.py}. The protocol repeats the $192^2\times201$, $240$-event experiment across three independent random seeds and reports aggregate gains over low-pass residual assimilation. At $3.0\%$ observation noise, OSNR reaches a mean trajectory PSNR of $67.0929$ dB versus $63.4397$ dB for low-pass, a mean gain of $3.6532$ dB with a worst-seed gain of $3.3038$ dB. At $5.0\%$ observation noise with the high-noise weak-atom smoothing profile, OSNR reaches $61.6653$ dB versus $59.6858$ dB for low-pass, a mean gain of $1.9795$ dB with a worst-seed gain of $1.7582$ dB. This establishes that the forced-flow advantage is not a single-seed artifact. \begin{table}[h] \centering \small \begin{tabular}{lccccc} \toprule Noise & Seeds & Mean recall & Mean OSNR & Mean low-pass & Worst gain \\ \midrule $3.0\%$ & $3$ & $73.33\%$ & $67.0929$ dB & $63.4397$ dB & $3.3038$ dB \\ $5.0\%$ & $3$ & $59.44\%$ & $61.6653$ dB & $59.6858$ dB & $1.7582$ dB \\ \bottomrule \end{tabular} \caption{Replicated high-noise forced Navier--Stokes sparse innovation assimilation at $192^2\times201$ with $240$ hidden forcing events. The reported gain is OSNR trajectory PSNR minus low-pass residual trajectory PSNR.} \label{tab:navier-stokes-high-noise-suite} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.78\linewidth]{../apps_industrial_breakthrough/navier_stokes_high_noise_suite_outputs_main/navier_stokes_high_noise_suite.png} \caption{Replicated high-noise forced-flow suite. OSNR remains ahead of smooth low-pass residual assimilation across all tested seeds at $3.0\%$ noise and after retuning the weak atom scale at $5.0\%$ noise.} \label{fig:navier-stokes-high-noise-suite} \end{figure} To probe whether the effect survives across a broader operating envelope, we added the ``destroyer'' matrix \texttt{apps\_industrial\_breakthrough/navier\_stokes\_destroyer\_protocol.py}. It evaluates $24$ forced-flow cases across four grid families ($96^2$, $128^2$, $160^2$, $192^2$), event counts from $48$ to $240$, noise levels from $1.0\%$ to $5.0\%$, and two random seeds per configuration. OSNR wins $22/24$ cases against the low-pass residual baseline, with mean trajectory gain $+4.2384$ dB and mean event recall $74.24\%$. At $3.0\%$ noise, OSNR wins all $12/12$ cases with mean gain $+2.9985$ dB. At larger grids ($128^2$, $160^2$, and $192^2$), OSNR wins every tested case; the only two losses occur in the smallest $96^2$ grid with the densest $96$-event, $5.0\%$ noise setting, where event overlap exceeds the available spatial resolution. \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Slice & Cases & Wins & Mean gain & Minimum gain \\ \midrule All destroyer cases & $24$ & $22$ & $+4.2384$ dB & $-0.7136$ dB \\ $1.0\%$ noise & $4$ & $4$ & $+13.6026$ dB & $+12.7441$ dB \\ $3.0\%$ noise & $12$ & $12$ & $+2.9985$ dB & $+0.7201$ dB \\ $5.0\%$ noise & $8$ & $6$ & $+1.4161$ dB & $-0.7136$ dB \\ $128^2$--$192^2$ grids & $16$ & $16$ & $+4.3012$ dB & $+1.7768$ dB \\ \bottomrule \end{tabular} \caption{Destroyer forced Navier--Stokes sparse assimilation matrix. Gains are OSNR trajectory PSNR minus low-pass residual trajectory PSNR.} \label{tab:navier-stokes-destroyer-protocol} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/navier_stokes_destroyer_protocol_outputs_main/navier_stokes_destroyer_protocol.png} \caption{Destroyer matrix summary. Bars show mean OSNR gain over low-pass residual forcing for each grid/event/noise configuration; labels show mean event recall.} \label{fig:navier-stokes-destroyer-protocol} \end{figure} \subsection{External PDEBench/FNO weather and fluid assimilation audits} \paragraph{Test~28 stabilizer boundary.} To move beyond internally generated forced-flow fields, we audited the hosted prediction tensors from the external \texttt{pdebench-fno-audit/fno-predictions} artifact. The target case is Test~28, a $512^2$ incompressible Navier--Stokes vorticity--Poisson benchmark. Each chunk stores FNO vorticity predictions $\hat\omega$, target vorticity $\omega$, and the published velocity-space nRMSE obtained by solving the Dirichlet Poisson problem \[ -\Delta\psi=\omega,\qquad \mathbf{v}=(\partial_y\psi,-\partial_x\psi), \] then comparing velocity fields after the first ten input frames. We implemented the same DST-I Poisson recovery in \texttt{apps\_industrial\_breakthrough/pdebench\_fno\_test28\_stabilizer.py} and verified that the recomputed FNO velocity nRMSE on chunk~00 matches the stored metric to within expected numerical drift. The OSNR diagnostic applies a deterministic spectral-viscosity operator to the FNO vorticity field, \[ \hat\omega_{\mathrm{osnr}} = \mathcal{F}^{-1}\!\left[M_K(k_x,k_y)\mathcal{F}\hat\omega\right], \] where $M_K$ is a compact rectangular frequency support. This is a blind post-processing stabilizer: it does not use held-out targets. We also tested a diagonal spectral transfer calibrated from five samples and applied to the remaining five samples; this is reported only as an assimilation diagnostic because it uses calibration targets. \begin{table}[h] \centering \small \begin{tabular}{lcc} \toprule Profile & Vorticity nRMSE & Velocity nRMSE \\ \midrule FNO artifact, chunk~00 & $1.390145$ & $0.243828$ \\ OSNR spectral viscosity, $K=64$ & $0.772164$ & $0.243647$ \\ OSNR best vorticity filter, $K=16$ & $0.542349$ & $0.244357$ \\ Diagonal spectral calibration, held-out & $0.497082$ & $0.266637$ \\ \midrule Oracle replace low modes, $K=8$ & -- & $0.057275$ \\ Oracle replace low modes, $K=32$ & -- & $0.023451$ \\ Oracle replace high modes, $K=32$ & -- & $0.243831$ \\ \bottomrule \end{tabular} \caption{External PDEBench/FNO Test-28 stabilizer audit on chunk~00. The OSNR spectral operator strongly suppresses vorticity outliers, but the official velocity-space metric changes only marginally and aggressive vorticity filtering can hurt velocity. The calibrated diagonal transfer is not a blind forecast result and is included to expose the metric boundary.} \label{tab:pdebench-fno-test28-stabilizer} \end{table} \begin{figure}[h] \centering \includegraphics[width=0.92\linewidth]{../apps_industrial_breakthrough/pdebench_fno_test28_outputs/test28_vorticity_panel.png} \caption{External PDEBench Test-28 vorticity panel for a held-out chunk sample. Spectral OSNR filtering removes large high-frequency FNO vorticity spikes, but this does not automatically translate into a large improvement in the benchmark velocity-space nRMSE.} \label{fig:pdebench-fno-test28-stabilizer} \end{figure} This external audit is a useful boundary result. It confirms that operator-spline spectral structure can repair raw vorticity instability in a real hosted FNO artifact, but it also prevents an overclaim: the official velocity metric is dominated by low-frequency phase and Poisson-integrated velocity structure. The oracle rows make this precise. Replacing only the lowest spectral modes of the FNO prediction with the target reduces velocity nRMSE from $0.243828$ to $0.057275$ at $K=8$ and $0.023451$ at $K=32$, whereas replacing high modes while leaving the low modes unchanged barely moves the metric. A SOTA-facing improvement on this benchmark therefore requires a velocity-aware low-mode dynamics corrector or Poisson-adjoint training objective, with sparse OSNR machinery reserved for high-frequency vorticity stabilization. We tested three follow-up low-mode correction families on held-out chunk~02 after calibrating on chunks~00--01. A diagonal vorticity-space transfer improved held-out vorticity nRMSE from $1.3716$ to $1.0479$ but worsened velocity nRMSE from $0.2303$ to $0.2491$. Direct velocity-space diagonal and mean-residual transfers also worsened the held-out metric, reaching best velocity nRMSE $0.2436$. Finally, a compact MPS-trained low-resolution velocity CNN fit the training loss but evaluated at $0.2488$ nRMSE on chunk~02. These negative results are informative: the low-mode error is not a stationary spectral bias and not solved by a small framewise image corrector. It is a sample-specific dynamical phase error. The next external benchmark attempt must either learn a genuine temporal low-mode evolution operator from the input history or select a benchmark where sparse/operator innovations, not phase drift, dominate the published metric. \paragraph{Sparse-station assimilation.} We therefore reframed the external FNO artifacts as sparse-station data assimilation problems, which is closer to operational weather and fluid monitoring. In \texttt{apps\_industrial\_breakthrough/pdebench\_weather\_sparse\_station\_assimilation.py}, the FNO forecast is treated as a neural dynamical prior. At each forecast step, a small number of station observations are used to solve either a closed-form per-channel affine correction or a low-rank DCT residual correction. On the three external Test~29 four-channel forecast configurations, this post-processing consistently improves the FNO forecast. The strongest case, \texttt{M01\_Eta01}, drops from mean nRMSE $0.00461$ to $0.000947$ with a rank-$8$ DCT residual fit from $512$ stations, a $79.5\%$ relative reduction. The \texttt{M10\_Eta01} case drops from $0.00551$ to $0.00159$ ($71.1\%$ reduction), while the harder \texttt{M10\_Eta001} case drops from $0.01204$ to $0.00908$ ($24.6\%$ reduction). We also tested the harder $512^2\times101$ vorticity artifact in \texttt{apps\_industrial\_breakthrough/pdebench\_vorticity\_sparse\_station\_assimilation.py}. Per-step affine station calibration reduces mean full-window vorticity nRMSE from $1.1802$ to $0.6872$ with $1024$ sparse stations, a $41.8\%$ relative reduction, but the velocity-space audit above shows that affine-only vorticity calibration is not the right final correction family. \begin{table}[h] \centering \small \begin{tabular}{lccc} \toprule External forecast artifact & FNO mean nRMSE & Sparse-station OSNR nRMSE & Relative gain \\ \midrule Test~29 \texttt{M01\_Eta01} & $0.004611$ & $0.000947$ & $79.5\%$ \\ Test~29 \texttt{M10\_Eta01} & $0.005507$ & $0.001594$ & $71.1\%$ \\ Test~29 \texttt{M10\_Eta001} & $0.012038$ & $0.009076$ & $24.6\%$ \\ Test~28 vorticity chunk~00 & $1.180217$ & $0.687232$ & $41.8\%$ \\ \bottomrule \end{tabular} \caption{External PDEBench/FNO sparse-station assimilation. The neural forecast is kept as the dynamical prior, while operator/dictionary corrections are solved from sparse observations without backpropagation.} \label{tab:pdebench-weather-station-assimilation} \end{table} The stronger Test~28 result comes from making the station correction operator-aware in the published metric. The runner \texttt{apps\_industrial\_breakthrough/pdebench\_vorticity\_dct\_station\_assimilation.py} fits, at each forecast step, an affine vorticity calibration followed by a rank-$32$ DCT residual from $1024$ contemporaneous vorticity stations. Unlike the affine-only station correction, this low-mode residual directly repairs the Poisson-integrated velocity structure. On all three available Test~28 chunks, the fixed profile improves both vorticity and the recomputed velocity metric: \begin{table}[h] \centering \small \begin{tabular}{lcccc} \toprule Chunk & FNO velocity & DCT-station velocity & FNO vorticity & DCT-station vorticity \\ \midrule 00 & $0.243826$ & $0.191702$ & $1.389743$ & $0.429227$ \\ 01 & $0.247869$ & $0.188312$ & $1.421356$ & $0.408249$ \\ 02 & $0.230292$ & $0.189971$ & $1.371190$ & $0.435880$ \\ \midrule Mean & $0.240662$ & $0.189995$ & $1.394096$ & $0.424452$ \\ \bottomrule \end{tabular} \caption{External PDEBench/FNO Test~28 DCT station assimilation. Metrics are computed after the first ten input frames. The correction uses $1024$ lattice stations, affine calibration, and a rank-$32$ DCT residual solve at each forecast step. Mean velocity nRMSE drops by $21.0\%$ and mean vorticity nRMSE drops by $69.5\%$ across the three local chunks.} \label{tab:pdebench-vorticity-dct-stations} \end{table} This per-frame result was the first Test~28 improvement in the project that moved the official velocity-space metric rather than only suppressing raw vorticity outliers. The next run added temporal structure to the station adapter. In \texttt{apps\_industrial\_breakthrough/pdebench\_vorticity\_temporal\_osnr\_station\_rescue.py}, chunk~00 selects the profile and chunks~01--02 are held out. For each trajectory, the frozen FNO rollout remains the neural dynamical prior, sparse contemporary vorticity stations are observed at each future frame, and one separable spatiotemporal OSNR residual is fitted over the whole forecast window, \[ \hat\omega(y,x,t)=a_t\hat\omega_{\mathrm{FNO}}(y,x,t)+b_t +\sum_{p,q,r} c_{pqr}\,\phi_p(y)\phi_q(x)\tau_r(t), \] where $\phi$ are spatial DCT atoms and $\tau$ is a temporal DCT basis over the post-input frames. The metric is again the Poisson-recovered velocity nRMSE plus raw vorticity nRMSE. \begin{table}[h] \centering \scriptsize \begin{tabular}{lccccc} \toprule Method on held-out chunks~01--02 & Stations/frame & Velocity nRMSE & Gain & Vorticity nRMSE & Gain \\ \midrule Frozen FNO prior & $0$ & $0.238631$ & -- & $1.396631$ & -- \\ Previous per-frame DCT, rank $32$ & $1024$ & $0.189142$ & $20.7\%$ & $0.422065$ & $69.8\%$ \\ Temporal OSNR, rank $32\times8$ & $2048$ & $0.112993$ & $52.6\%$ & $0.349838$ & $75.0\%$ \\ Temporal OSNR, rank $32\times12$ & $4096$ & $0.106656$ & $55.3\%$ & $0.328363$ & $76.5\%$ \\ Temporal OSNR, rank $48\times12$ & $4096$ & $0.103622$ & $56.6\%$ & $0.327553$ & $76.5\%$ \\ \bottomrule \end{tabular} \caption{External PDEBench/FNO Test~28 temporal OSNR station rescue. The final row uses $4096/262144=1.5625\%$ of grid sites per frame and fits one spatiotemporal residual over each forecast trajectory. Hyperparameters are selected on chunk~00 and reported on held-out chunks~01--02.} \label{tab:pdebench-vorticity-temporal-stations} \end{table} This became the strongest closed-form Test~28 external FNO result in the workspace. Across all three local chunks, the rank-$48\times12$ temporal adapter reduces mean velocity nRMSE from $0.240406$ to $0.103368$ and mean vorticity nRMSE from $1.394469$ to $0.330206$. It is still an assimilation result, not a blind forecast: the method uses contemporary sparse measurements. That distinction is important, but it is also exactly the operational setting where station, buoy, radar, and satellite observations are available and a neural forecast acts as the dynamical prior. The broader mechanism is now clearer than in the first Test~28 audit: OSNR can act as a closed-form, low-rank test-time correction layer on top of a frozen neural PDE forecaster, and the same plug-in idea already transferred from external Darcy sparse assimilation to time-dependent vorticity forecasts. We then ran the harder SOTA-facing comparison: a trained sparse neural assimilator under the same station protocol. The runner \texttt{apps\_industrial\_breakthrough/pdebench\_vorticity\_neural\_sparse\_assimilation\_baseline.py} trains on chunks~00--01 and reports held-out chunk~02. Its inputs are the frozen FNO vorticity frame, the same-station affine calibration, sparse target and residual maps, a station mask, coordinates, and forecast time. A compact $770{,}241$-parameter U-Net predicts a $128^2$ residual that is upsampled to the full $512^2$ grid before both vorticity and Poisson velocity metrics are computed. This is not a no-backprop OSNR result; it is the competent trained sparse neural baseline that the closed-form adapter must be compared against. \begin{table}[h] \centering \scriptsize \begin{tabular}{lcccc} \toprule Method on held-out chunk~02 & Stations/frame & Velocity nRMSE & Gain vs FNO & Vorticity nRMSE \\ \midrule Frozen FNO prior & $0$ & $0.230295$ & -- & $1.371558$ \\ FNO + temporal OSNR & $256$ & $0.198508$ & $13.8\%$ & $0.551600$ \\ Sparse neural assimilator & $256$ & $0.073335$ & $68.2\%$ & $0.191289$ \\ FNO + temporal OSNR & $512$ & $0.147680$ & $35.9\%$ & $0.461464$ \\ Sparse neural assimilator & $512$ & $0.031716$ & $86.2\%$ & $0.154792$ \\ FNO + temporal OSNR & $1024$ & $0.134717$ & $41.5\%$ & $0.416209$ \\ Sparse neural assimilator & $1024$ & $0.019973$ & $91.3\%$ & $0.139996$ \\ Sparse neural assimilator, repeat seed & $1024$ & $0.021120$ & $90.8\%$ & $0.140282$ \\ FNO + temporal OSNR & $4096$ & $0.110239$ & $52.1\%$ & $0.338437$ \\ Sparse neural assimilator & $4096$ & $\mathbf{0.015279}$ & $\mathbf{93.4\%}$ & $0.155159$ \\ \bottomrule \end{tabular} \caption{PDEBench Test~28 trained sparse neural assimilation on held-out chunk~02. The $1024$-station row uses only $1024/262144=0.390625\%$ of grid sites per frame and is stable under one seed repeat. The neural model is trained with backpropagation and is included as the relevant SOTA-facing sparse-assimilation comparator.} \label{tab:pdebench-vorticity-neural-stations} \end{table} This shifts the Test~28 frontier. The $1024$-station neural row reduces velocity error by about $91\%$ against the frozen FNO and by about $85\%$ against the matched closed-form FNO+OSNR row. The $4096$-station row reaches the best velocity value, $0.015279$, but the $1024$ row is the cleaner observation-efficiency result. The fixed FNO-tuned temporal OSNR adapter does not transfer unchanged onto the trained neural prior: at $1024$ stations it worsens the neural velocity row from $0.019973$ to $0.075423$, and at $4096$ from $0.015279$ to $0.048526$. This negative adapter result is useful. Once the neural model has learned the low-frequency station-conditioned correction, the next OSNR layer must be selected specifically for a neural prior, with identity/gating/high-ridge/low-rank candidates or an orthogonalized residual space. Reusing the FNO prior's adapter is not valid. The active-station follow-up then exposed that the original closed-form gap was mostly geometric. We extended the same runner with centered, interior, space-filling, and rounded phase-lattice station policies while keeping the train/test split, neural architecture, epochs, and Poisson velocity metric fixed. The old lattice includes boundary-heavy samples; a centered or interiorized lattice spends the same budget on Fourier-compatible interior coverage. Table~\ref{tab:pdebench-vorticity-active-stations} shows the result on held-out chunk~02. \begin{table}[h] \centering \scriptsize \begin{tabular}{lccccc} \toprule Station policy & Stations & FNO+OSNR velocity & FNO+OSNR vorticity & Neural velocity & Neural vorticity \\ \midrule Old edge lattice & $512$ & $0.147680$ & $0.461464$ & $0.031716$ & $0.154792$ \\ Centered lattice & $512$ & $0.036161$ & $1.289455$ & $0.042066$ & $1.085325$ \\ Space filling & $512$ & $0.142717$ & $0.461791$ & $0.031003$ & $0.176815$ \\ \midrule Old edge lattice & $1024$ & $0.134717$ & $0.416209$ & $0.019973$ & $0.139996$ \\ Centered lattice & $1024$ & $0.020786$ & $1.281390$ & $0.034668$ & $1.077794$ \\ Best rounded phase $(0.50,0.25)$ & $1024$ & $\mathbf{0.020563}$ & $1.285525$ & $0.038600$ & $1.078815$ \\ Space filling & $1024$ & $0.118799$ & $0.350398$ & $0.017219$ & $0.176978$ \\ \midrule Old edge lattice & $2048$ & $0.119788$ & $0.364925$ & $0.020052$ & $0.147939$ \\ Space filling & $2048$ & $0.108844$ & $0.336009$ & $0.016828$ & $0.212521$ \\ Old edge lattice & $4096$ & $0.110239$ & $0.338437$ & $\mathbf{0.015279}$ & $0.155159$ \\ Space filling & $4096$ & $0.106388$ & $0.340152$ & $0.015478$ & $0.221071$ \\ \bottomrule \end{tabular} \caption{PDEBench Test~28 active station-geometry audit on held-out chunk~02. Regular interior station geometry nearly closes the velocity gap between closed-form temporal OSNR and the trained sparse neural assimilator at $1024$ stations, but does not repair raw vorticity. Space-filling stations improve the trained neural velocity curve at $1024$--$2048$ stations but worsen vorticity and do not beat the old $4096$-station neural velocity frontier.} \label{tab:pdebench-vorticity-active-stations} \end{table} The best closed-form phase row reduces FNO velocity nRMSE from $0.230295$ to $0.020563$ using only $1024/262144=0.390625\%$ contemporary station sites per frame. This is an $84.7\%$ reduction relative to the old same-budget edge-lattice OSNR row and is only about $3\%$ worse than the trained $1024$-station neural row. It also beats one repeat seed of that neural row ($0.021120$). The caveat is just as important: the same regular/phase lattice rows leave vorticity near $1.28$, so the win is a Poisson-velocity low-mode correction, not full vorticity reconstruction. Space-filling gives the complementary behavior: the neural model reaches new $1024$- and $2048$-station velocity-efficiency rows, $0.017219$ and $0.016828$, but with worse vorticity and no improvement over the old $4096$-station velocity frontier. A local phase refinement around $(0.50,0.25)$ was locally saturated and quantized by integer-grid rounding. A naive one-mask hybrid that concatenates $50\%$ or $75\%$ centered/interior lattice stations with space-filling fill points was decisively negative: the best $1024$-station hybrid neural velocity was only $0.131524$, and closed-form hybrid velocity was worse than the frozen FNO. Thus the next Test~28 step should be a vorticity-aware two-geometry or two-head adapter: keep Fourier-compatible interior stations for the low-mode velocity correction, add a separate high-frequency/vorticity residual mechanism, and gate any OSNR residual against the trained neural prior rather than reusing the FNO-prior adapter blindly. The two-head follow-up made this decomposition explicit. Because the benchmark velocity is recovered by a DST-I Poisson solve, the useful fusion basis is not a generic DCT split but the same sine basis that diagonalizes the reported metric. The runner \texttt{pdebench\_vorticity\_two\_head\_frequency\_adapter.py} keeps the phase-lattice OSNR head for low modes, adds a separate space-filling OSNR vorticity head for high modes, selects the DST cutoff and scalar weight on chunk~01, and reports chunk~02. Table~\ref{tab:pdebench-vorticity-dst-two-head} summarizes the resulting Pareto rows. \begin{table}[h] \centering \scriptsize \begin{tabular}{lccc} \toprule Method & Station observations/frame & Velocity nRMSE & Vorticity nRMSE \\ \midrule FNO & $0$ & $0.230295$ & $1.371558$ \\ Phase low head & $1024$ & $0.020563$ & $1.285525$ \\ Space-filling high head & $1024$ & $0.118799$ & $0.350398$ \\ DST two-head & $1024+1024$ & $0.023612$ & $0.192818$ \\ DST two-head, larger high head & $1024+2048$ & $0.023712$ & $0.195319$ \\ Cached neural comparator & $1024$ & $0.017219$ & $0.176978$ \\ Corrected neural repeat & $1024$ & $0.017219$ & $0.176978$ \\ Neural + DST high-pass & $1024$ & $0.017060$ & $0.113794$ \\ Corrected neural repeat & $2048$ & $0.016828$ & $0.212521$ \\ Neural + DST high-pass & $2048$ & $0.016403$ & $0.097739$ \\ Corrected neural repeat & $4096$ & $0.015478$ & $0.221071$ \\ Neural + DST high-pass & $4096$ & $\mathbf{0.014927}$ & $\mathbf{0.097201}$ \\ Corrected neural repeat, seed $20260611$ & $1024$ & $0.020466$ & $0.177627$ \\ Neural + DST high-pass, seed $20260611$ & $1024$ & $0.020210$ & $0.095205$ \\ \bottomrule \end{tabular} \caption{PDEBench Test~28 DST two-head and neural-prior high-pass adapters on held-out chunk~02. Separate station counts indicate distinct low-mode and high-mode observation sets for the closed-form rows. Each neural high-pass row uses the same deterministic space-filling station set as its neural comparator, selects the DST split on chunk~01, and improves both reported metrics on chunk~02.} \label{tab:pdebench-vorticity-dst-two-head} \end{table} The closed-form DST row is the first Test~28 adapter in the workspace that substantially improves both sides of the earlier closed-form tradeoff: vorticity drops from the standalone space-filling value $0.350398$ to $0.192818$, while velocity remains close to the phase-lattice value ($0.023612$ versus $0.020563$). Increasing the high-frequency head to $2048$ stations does not improve the fused row, so the bottleneck is not simply high-head station count. The neural-prior row is the clean frontier. The first repeat used the correct station geometry but the wrong statistics-sampling seed; after matching the baseline convention, the runner exactly reproduces the cached $1024$-station neural row. A validation-selected DST high-pass OSNR layer with cutoff $80$ and weight $0.2$ then improves held-out velocity from $0.017219$ to $0.017060$ and raw vorticity from $0.176978$ to $0.113794$, using the same deterministic space-filling station set. The follow-up ladder strengthens the claim. At $2048$ stations, the vorticity-aware selector chooses cutoff $128$ and weight $0.1$, improving the neural row from $0.016828$/$0.212521$ to $0.016403$/$0.097739$. At $4096$ stations, cutoff $128$ and weight $0.05$ improve $0.015478$/$0.221071$ to $0.014927$/$0.097201$. A second $1024$-station seed repeats the pattern, improving $0.020466$/$0.177627$ to $0.020210$/$0.095205$. This reverses the earlier negative neural+OSNR result, where an FNO-tuned residual damaged the trained neural prior: the useful adapter is neural-prior-specific and orthogonalized into the DST high-frequency space. \begin{table}[h] \centering \scriptsize \begin{tabular}{lcccc} \toprule Gate row & Fit/gate stations & Choices & Neural vel/vort & Gated vel/vort \\ \midrule $1024$, seed $20260610$ & $768/256$ & $112{:}10,\ 128{:}0$ & $0.017219/0.176978$ & $0.016922/0.096622$ \\ $2048$, seed $20260610$ & $1536/512$ & $112{:}10,\ 128{:}0$ & $0.016828/0.212521$ & $0.016403/0.099564$ \\ $4096$, seed $20260610$ & $3072/1024$ & $112{:}6,\ 128{:}4$ & $0.015478/0.221071$ & $0.014927/0.098693$ \\ $1024$, seed $20260611$ & $768/256$ & $112{:}10,\ 128{:}0$ & $0.020466/0.177627$ & $0.020210/0.097617$ \\ \bottomrule \end{tabular} \caption{No-leakage station-heldout gate for the Test~28 neural-prior DST adapter. The gate chooses identity, cutoff $112$, or cutoff $128$ from held-out station residuals only, then refits the chosen correction on all available stations for the reported field. Identity is never selected in these runs.} \label{tab:pdebench-vorticity-station-gate} \end{table} The gate table removes the remaining hand-picked-selector weakness. The decision uses only contemporary station values: $25\%$ of stations are withheld from the gate fit, candidate residuals are scored on those stations, and the selected candidate is then refit on the full station set. This standard cross-validation pattern is operationally different from using full-field validation metrics. The gate is conservative relative to the full-field oracle, which would choose cutoff $128$ in all four runs; it often chooses cutoff $112$ instead. The price is a small vorticity gap versus the oracle, but the no-leakage rows still improve both velocity and vorticity over the trained neural prior at every tested budget and on the second seed. A fit-only ablation that permanently discards the gate stations damages the velocity metric, so the deployable protocol is station-heldout selection followed by all-station refit. The more aggressive follow-up removes the assimilation advantage entirely. The blind diffusion-refiner runner sees only the first ten true vorticity frames of each held-out Test~28 trajectory at inference time. It never reads the FNO prediction, never observes future stations, and uses the cached FNO tensor only as an evaluation comparator. The model treats forecasting as an iterative refinement-time PDE, \[ u_{k+1}=P\!\left(u_k+\eta F_\theta(u_k,\mathrm{history},\tau,k)\right), \] where $P$ is either the identity or a DST spectral-viscosity projection. After the first run showed a clean failure mode--excellent vorticity but weaker Poisson velocity--we added a light Poisson-weighted DST coefficient loss and then a differentiable low-mode Poisson-velocity loss. \begin{table}[h] \centering \scriptsize \begin{tabular}{lcc} \toprule Blind-from-history row on chunk~02 & Velocity nRMSE & Vorticity nRMSE \\ \midrule FNO comparator & $0.230295$ & $1.371558$ \\ Persistence from frame 9 & $0.693986$ & $1.175046$ \\ Linear extrapolation from frames 8/9 & $1.497893$ & $4.620379$ \\ Pure learned refiner, $r128/e10$ & $0.303547$ & $0.414996$ \\ DST projected refiner, $r128/e10$ & $0.297367$ & $0.415566$ \\ DST + Poisson loss, seed $20260612$ & $0.261985$ & $0.405375$ \\ DST + Poisson loss, seed $20260613$ & $0.261936$ & $0.400889$ \\ DST + Poisson loss, seed $20260614$ & $0.252706$ & $0.396079$ \\ DST + Poisson loss, seed $20260615$ & $0.281126$ & $0.409406$ \\ DST + Poisson + velocity loss, seed $20260614$ & $0.246451$ & $0.396002$ \\ DST + Poisson + velocity loss, seed $20260613$ & $0.265005$ & $0.405434$ \\ DST + Poisson + velocity loss, seed $20260612$ & $0.268111$ & $0.413803$ \\ Uniform ensemble, seeds $20260614/13/12$ & $0.233286$ & $0.374172$ \\ Uniform ensemble, seeds $20260614/13/12/11/15$ & $\mathbf{0.228610}$ & $\mathbf{0.366712}$ \\ Validation-locked subset, seeds $20260614/11/15$ & $0.233002$ & $0.370837$ \\ Validation-locked weighted, seeds $20260614/12/11/15$ & $0.231505$ & $0.370884$ \\ Validation-locked weighted, seeds $20260614/13/12/11/15$ & $0.231374$ & $0.371718$ \\ \bottomrule \end{tabular} \caption{Blind Test~28 from-initial-history refiner. The OSNR/DST rows do not use FNO predictions or future observations at inference time. They are trained from chunks~00/01 and evaluated on held-out chunk~02; the FNO row is a frozen external comparator.} \label{tab:pdebench-vorticity-blind-refiner} \end{table} The direct velocity objective tightened the single-seed frontier from $0.252706$ to $0.246451$ while preserving the vorticity win. More importantly, seed diversity exposed an ensemble effect rather than a single lucky run. A three-seed uniform average nearly closes the FNO velocity gap, and the five-seed uniform ensemble becomes the first blind from-initial-history Test~28 row in this project to beat the hosted FNO comparator on both reported metrics: velocity improves from $0.230295$ to $0.228610$ ($0.73\%$), while vorticity drops from $1.371558$ to $0.366712$ ($73.3\%$). We then froze an explicit validation-locked model-selection protocol in \texttt{pdebench\_vorticity\_blind\_locked\_ensemble.py}: train candidate seeds on chunk~00, select a uniform seed subset on validation chunk~01 with the predeclared score velocity plus $0.02$ times vorticity, refit only the selected seeds on chunks~00/01, and evaluate chunk~02 once. That protocol selects seeds $20260614/11/15$ and reaches $0.233002$ velocity and $0.370837$ vorticity on held-out chunk~02. The locked row therefore confirms the blind raw-vorticity result--$73.0\%$ lower vorticity than FNO--but it does not yet confirm the exploratory velocity edge, missing FNO velocity by $1.18\%$. The scientific signal is sharp but narrower than the frontier row: without FNO input or future observations, OSNR/DST refinement has a defensible no-leakage mechanism for collapsing raw vorticity, while the Poisson-velocity win still requires a better locked selector or dynamics model. We then tested whether the missing velocity margin was simply a selector issue. The weighted locked runner, \texttt{pdebench\_vorticity\_blind\_locked\_weighted\_ensemble.py}, builds validation-only low-resolution prediction quadratics and searches a convex $0.05$ simplex grid with the predeclared score mean velocity plus $0.02$ times mean vorticity plus $0.25$ times velocity p90 plus $0.05$ times maximum velocity. Allowing four or five active seeds selects weights $(0.20,0.15,0.20,0.45)$ on seeds $20260614/12/11/15$ and reaches $0.231505$ velocity, $0.370884$ vorticity on chunk~02. Forcing all five seeds active selects weights $(0.20,0.05,0.10,0.15,0.50)$ and reaches the best strict locked velocity row so far: $0.231374$ velocity and $0.371718$ vorticity. This narrows the locked velocity gap from $1.18\%$ to $0.47\%$ relative to FNO, but still does not beat the FNO velocity comparator. The vorticity result remains stable at roughly $72.9\%$ lower error. Thus the next velocity gain is unlikely to come from selector-only sweeps; it likely requires a stronger cross-mode or horizon-conditioned dynamics model. The negative controls are also informative. A no-backprop per-mode DST ridge dynamics runner reaches only $0.352961$ velocity and $0.555332$ vorticity at $r128/K64$; the validation-selected cross-mode kernel dynamics follow-up is worse still at $0.868662$ velocity and $0.887061$ vorticity; naive $256$-grid scaling reaches only $0.320404$ velocity and $0.462010$ vorticity after rescaling the Poisson loss; and two $20$-epoch schedules improve training loss without improving held-out velocity ($0.270341$ and $0.267655$). Thus the current blind lesson is specific: iterative learned refinement plus light Poisson-weighted OSNR/DST structure and seed diversity are the useful branch, while per-mode closed-form spectral extrapolation, naive cross-mode kernels, naive resolution scaling, and longer single-seed training do not close the single-seed velocity gap. This is not a claim of beating end-to-end global weather systems such as GraphCast or GenCast. It is a more precise and immediately defensible claim: operator-spline station assimilation can dramatically improve external neural PDE forecasts when a small number of contemporary observations are available. That setting is scientifically meaningful because real forecasting systems already assimilate sparse stations, buoys, sondes, radar, and satellite products; the OSNR contribution is a closed-form, low-rank correction layer that can sit on top of a neural forecaster without retraining it. We also tested a local compact-kernel variant in \texttt{apps\_industrial\_breakthrough/pdebench\_weather\_rbf\_station\_assimilation.py}. Here the station residual after affine calibration is interpolated through a Gaussian RBF kernel, which better matches spatially local forecast-error structure than a global DCT basis at low station counts. On \texttt{M01\_Eta01}, affine+RBF improves the best result further from $0.000947$ to $0.000819$, an $82.2\%$ reduction from the FNO baseline. On \texttt{M10\_Eta01}, RBF reaches $0.001640$ and already obtains a $52.9\%$ reduction with only $16$ stations, while the global DCT correction remains slightly better at the highest station budget. On \texttt{M10\_Eta001}, DCT remains the better choice. The conclusion is algorithmic rather than cosmetic: the assimilation layer should select its correction dictionary from the forecast-error geometry. Smooth global biases favor low-rank DCT; localized residual structure favors compact RBF/operator-spline kernels. The follow-up hybrid runner \texttt{apps\_industrial\_breakthrough/pdebench\_weather\_hybrid\_station\_assimilation.py} combines the two correction families in one sequential closed-form layer: affine calibration, rank-$8$ DCT residual fitting, and station-centered Gaussian RBF residual interpolation. A low-budget targeted sweep shows that blindly mixing dictionaries can underperform the best single family because station equations are split across redundant atoms. At the higher $512$-station random budget, however, the hybrid layer improves all three external forecasts: \texttt{M01\_Eta01} drops to $0.000707$ ($84.7\%$ reduction), \texttt{M10\_Eta01} drops to $0.001438$ ($73.9\%$ reduction), and \texttt{M10\_Eta001} drops to $0.008584$ ($28.7\%$ reduction). This became the random-station reference for the adaptive placement tests. The design rule is now clearer: use a single compact dictionary at very sparse station counts, but switch to a hybrid global--local operator dictionary once the observation budget is high enough to identify both smooth bias and localized forecast residuals. \paragraph{Adaptive Test~29 station placement.} The next experiment, \texttt{apps\_industrial\_breakthrough/pdebench\_weather\_adaptive\_station\_assimilation.py}, tests whether station placement can reduce the observation budget. The correction dictionary is kept fixed and only the station policy changes. Forecast-gradient, DCT-leverage, and two-stage pilot-residual policies are not reliable: they oversample high-variation or high-error regions and leave the RBF/DCT normal equations poorly covered. A farthest-point space-filling policy is better because it improves interpolation coverage and dictionary conditioning. The first adaptive sweep nearly matched the $512$-random hybrid with $256$ stations. A targeted follow-up then retuned only the RBF length scale and showed that $256$ space-filling stations are enough to beat the $512$-random reference on all three external Test~29 files: \texttt{M01\_Eta01} reaches $0.000705$, \texttt{M10\_Eta01} reaches $0.001377$, and \texttt{M10\_Eta001} reaches $0.007649$. The same policy at $288$ stations improves all three again, and $384$ stations gives the strongest current Test~29 assimilation results: $0.000646$, $0.001267$, and $0.007274$, respectively. Thus the useful high-information criterion for this external artifact is not local forecast-error magnitude; it is conditioned spatial coverage for the global--local operator dictionary. The result is a concrete observation-efficiency gain: half as many stations now beat the former $512$-station hybrid reference, without retraining the FNO and without backpropagating through the assimilation layer. For reproducibility, this weather result should be read as a frozen-neural-prior plus deterministic assimilation architecture, not as a newly trained neural network. The neural component is the external FNO prediction tensor already present in the hosted audit artifact. The OSNR code never changes the FNO weights and does not run reverse-mode differentiation. Each Test~29 file contains arrays \texttt{preds} and \texttt{targets} with shape $10\times128\times128\times21\times4$: ten held-out forecast samples, a $128^2$ spatial grid, $21$ forecast frames, and four weather/PDE channels. For each sample $n$, time index $\tau$, and channel $c$, the layer treats the FNO forecast $p_c(x)$ as a dynamical prior and receives contemporary station observations $y_c(x_i)$ at a selected set $S$ of grid sites. The correction is the three-block operator \[ \mathcal{A}_{S,r,\ell}(p,y) = \mathcal{R}_{S,\ell}\!\left( \mathcal{D}_{S,r}\!\left( \mathcal{C}_{S}(p,y),y\right),y\right), \] where $\mathcal{C}_S$ is per-channel affine calibration, $\mathcal{D}_{S,r}$ is a low-rank DCT residual solve, and $\mathcal{R}_{S,\ell}$ is a compact Gaussian-RBF station residual solve. The first block solves \[ (a_c,b_c)=\arg\min_{a,b}\sum_{i\in S}(a\,p_c(x_i)+b-y_c(x_i))^2, \qquad q_c^{(0)}(x)=a_c p_c(x)+b_c. \] The DCT block builds a tensor-product dictionary $\Phi_r\in\mathbb{R}^{HW\times r^2}$ with normalized atoms \[ \phi_{k_y,k_x}(i,j) = Z^{-1}_{k_y,k_x} \cos\!\left(\frac{\pi(i+1/2)k_y}{H}\right) \cos\!\left(\frac{\pi(j+1/2)k_x}{W}\right), \qquad 0\leq k_y,k_x