Variable-permittivity edge propagator, and a correction to how #52 was framed - #55
Open
arnaudon wants to merge 28 commits into
Open
Variable-permittivity edge propagator, and a correction to how #52 was framed#55arnaudon wants to merge 28 commits into
arnaudon wants to merge 28 commits into
Conversation
Groundwork for the within-edge resolution problem, and a correction to how #52 framed it. I had recommended computing the hole-burning integral analytically per edge, "as compute_mode_competition_matrix already does". That is wrong: there is no such integral in the operator. construct_weight_matrix uses one k_e per edge, and 1/(exp(2i k_e l_e) - 1) is the exact solution *only* for constant eps. Spatial hole burning makes eps vary within an edge, which breaks the piecewise-constant assumption the quantum-graph method is built on. Oversampling restores it by subdividing, and pays in the size of the eigenproblem -- 76803 nodes against 243 edges on the buffon. An edge with varying eps still has an exact 2x2 transfer matrix. Computing it per edge keeps the eigenproblem at its original size and makes the sub-interval count local, independent and parallel. Two measured results, in examples/audit/probe_edge_propagator.py against a DOP853 reference at rtol 1e-13 on a buffon-like edge (~28 oscillations, 4% saturation ripple): * Oversampling *is* second-order Magnus. Freezing eps at each sub-edge midpoint and multiplying constant-eps propagators is exactly exp(h A(x_mid)) for this system; they agree to round-off. So the current scheme is O(h^2) and no budget tuning improves that. * Fourth-order Magnus needs 6400 sub-intervals for 1e-8 where second order needs 819200 -- 128x fewer, on top of taking the count out of the matrix. The propagator reproduces the closed form exactly for constant eps, so it can replace it without changing passive behaviour; that is asserted in the tests along with the O(h^2) vs O(h^4) orders, Wronskian conservation, and the magnus2-equals-piecewise-constant identity. Wiring it into construct_incidence_matrix / construct_weight_matrix -- where the open and directed boundary models need handling -- is the follow-up. 190 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Second step toward #52/#53. The secular matrix L = B^T W^-1 B turns out to be exactly a sum of independent per-edge 2x2 Dirichlet-to-Neumann blocks, and the DtN map of a varying-eps edge is available from its transfer matrix. So the matrix can stay one-node-per-vertex while the within-edge resolution moves into local per-edge work -- which is the whole reason oversampling is unaffordable on the buffon (76803 nodes against 243 edges). Derivation, with M the edge transfer matrix and det M = 1 (Wronskian): psi'_u = (psi_v - M11 psi_u)/M12, M21 - M22 M11/M12 = -1/M12 => D = (1/M12) [[M11, -1], [-1, M22]], L_edge = -i D For constant eps this collapses to -i q cot(ql) on the diagonal and i q/sin(ql) off it, i.e. exactly what construct_laplacian already builds -- asserted to 1e-12 on both the closed and open models. That exact reduction is the property that makes it safe to swap in. Boundary edges under open_model="open" carry the outgoing-wave block q/(e^2-1) [[1, -e], [-e, 1]] instead, which is only valid for constant eps. Those edges are the passive leads (pump = 0), so hole burning leaves them uniform; the function keeps the closed form for them and *raises* if handed a varying profile there, rather than returning a matrix that is not the DtN map of anything. The directed models zero a different set of B/B^T entries and are rejected explicitly for the same reason. Not yet wired into the solver -- that is the next step, and what actually closes the issues. 198 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Third step toward #52/#53. The saturated permittivity is built from |E(x)|^2 *within* the edge, so the solver needs the field there, not just at the nodes. edge_step_propagators exposes the per-sub-interval factors (whose ordered product is the transfer matrix), and edge_field_samples propagates (psi, psi') from one end to give psi on the sub-interval grid. Keeping the factors rather than re-propagating avoids doing the work twice: the same sweep that builds the edge's DtN block also yields the field the next profile is built from. Verified against the analytic cos(qx) for constant eps, and against the transfer matrix at the endpoints for a varying one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The load-bearing check for #52/#53. On a 5-edge ring with a 10% standing-wave ripple in eps, the oversampled answer converges to the per-edge-DtN answer at a clean O(h^2) -- 6.24e-4, 1.61e-4, 4.04e-5 at 16/32/64 sub-edges per edge, ratios 3.88 and 3.99 -- so the two are solving the same problem. The DtN operator reaches it directly at |lambda_min| = 1.1e-13 with a 5x5 matrix against 320x320. The ripple shifts the mode by 1.8e-1, four orders above the agreement, so the comparison is not vacuous. Two traps found while writing this, both of which produced convincing nonsense before being caught, and both now guarded: * A minimum on the bracket boundary is not a mode. The first version took a fixed window around a guess; every configuration returned its own endpoint and they agreed to 1e-13 with none having found a root. find_mode now requires an interior minimum and raises otherwise. * Sub-edge eps must be placed geometrically. Accumulating sub-edge lengths in work.edges order assumes that order walks each parent edge end to end; it does not, and the scrambled profile put the reference on a different mode (6.42 against 5.42). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…ion it needs Fourth step toward #52/#53: the hole-burning physics carried on the varying operator instead of an oversampled graph. node_solution_varying gives the null vector at real k; edge_field_profiles propagates (psi, psi') from each edge's start -- taking psi'_u = (psi_v - M11 psi_u)/M12 from the same transfer matrix that built the operator's block -- to sample |E(x)|^2 *inside* the edge; saturated_eps_profiles turns those into per-edge eps(x); salt_residuals_varying is the acceptance test. With the saturation switched off (a = 0) this agrees with modes.salt_residuals to 1.5e-15, which exercises the whole chain: operator, gain, boundary handling, profile construction. Two things had to be got right, both found by measurement rather than reasoning: * The pump norm is *unconjugated*. modes._graph_norm contracts with node_solution.T (not .conj().T) against a compute_z_matrix built from exp(2i l k)/(2i k), i.e. it integrates psi^2, while the per-edge intensity it divides is |E|^2 (_mean_intensity_from_flux uses k + conj(k)). That asymmetry is the SALT biorthogonal convention. Using int|psi|^2 instead left a constant 2.2% offset on every edge -- invisible at a = 0, and a silent rescaling of what an amplitude means at a > 0. Matching it brings the per-edge means to within 1e-6 of the existing path. * The norm quadrature must be Simpson, not trapezoid. The propagator is fourth-order, so an O(h^2) quadrature would have set the accuracy of the whole profile and thrown away what the Magnus step bought (means converging 4.4e-5 -> 1.1e-5 -> 2.8e-6 over doublings). Not yet the solver loop, so the issues are still open. Note also that comparing the two paths at an arbitrary (D0, a) is *not* a valid check: such a state is not a SALT solution, |lambda_1| never approaches zero, and the argmin of a non-vanishing residual is not a physical quantity. That is recorded in the test class so the next person does not repeat it. 209 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…rsampled Fifth step toward #52/#53. solve_salt_varying is the varying-operator counterpart of modes.solve_salt_fixed_set: per mode the unknowns are (k, a) and the equations are Re lambda_1(k) = Im lambda_1(k) = 0 -- square -- with the hole-burning fields frozen during each least-squares solve and refreshed after. It converges to genuine solutions on an 8-node Fabry-Perot that is never oversampled: D0/D0_thr k a residual iters 1.2 11.41362806 0.08610 5.5e-07 10 1.5 10.44115645 0.33935 7.8e-07 10 2.0 9.46923150 0.54464 9.9e-07 11 Two things worth recording, both mistakes in my own testing rather than the code: * The first threshold scan ran D0 only to 0.12 and concluded "no threshold". |lambda_1| was falling linearly (3.18 -> 2.38) and crosses zero near 0.41; the scan simply stopped short. The amplitudes correctly collapsed to zero throughout, which is what below-threshold should do. * The k values above jump between modes because the starting guesses came from a scan that wanders. Each solve converges to a genuine root; which root depends on where it starts. Pump continuation is the fix, and is what the existing _full_salt_newton_impl already does. With D0 = 0 the operator has no gain, so the residual has no gradient in `a` and least-squares leaves it at its initial value. That is correct for an unconstrained sub-problem, and the result is flagged rather than dressed up: converged is False and the residual stays large. Pinned in a test, since the tempting assertion (a collapses to 0) is wrong. 211 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The direct answer to #53's complaint. Where the oversampled solver gave a non-monotone summed output (-40.8% across one pump step on mini_buffon) and residuals of 1e-2 deep above threshold, the same sweep through solve_salt_varying -- matrix fixed at 8x8, resolution in per-edge transfer matrices -- gives: 15 of 15 pumps converged, residuals 3e-7..9e-7, amplitude strictly monotone 0.037 -> 1.372, k stable to 1e-5 from 1.05x to 3.0x threshold, 11.4 s. The continuation is load-bearing, and the README says so: starting each pump from a scan-derived guess instead of the previous solution makes the solver land on different modes (k jumping 11.4 -> 10.4 -> 9.5), each a genuine root but not the same branch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The measurement that answers #52. Oversampling needs ~76600 nodes for lambda/12 on this graph and the default cap gives 0.47 samples per wavelength -- four times below Nyquist, so the within-edge field is aliased rather than coarse. On the per-edge-DtN operator the matrix stays 208x208 at every resolution: n_steps samples/wavelength matrix build |lambda| 64 2.3 208x208 0.29s 1.73132837 256 9.2 208x208 0.99s 1.73129737 d=1.29e-05 512 18.5 208x208 1.95s 1.73129710 d=6.46e-07 1024 36.9 208x208 3.96s 1.73129709 d=3.34e-08 18.5 samples per wavelength against 0.47, converged to 6e-7, in a matrix 368x smaller than the equivalent subdivision. The first version of this probe ran at zero amplitude, which leaves the saturation denominator at 1: the profile is then constant, the propagator exact at any n_steps, and |lambda| identical to 1e-12 across the sweep. It looked like flawless convergence while exercising none of the varying machinery. The script now uses a genuine saturated field and says so. AUDIT.md gains a section 9 recording that section 8's diagnosis -- a flat node budget -- was wrong, and that the real obstacle is the operator's piecewise-constant assumption. #52 and #53 stay open: the new path is public API but is not yet what intensity_method: full_salt_newton runs, and it has not been cross-checked against examples/audit/independent_salt/. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…arying
compute_modal_intensities_varying does pump continuation through
solve_salt_varying and returns the usual ("modal_intensities", D0) columns, so
`python -m netsalt lasing config.yaml` runs it. Verified end to end on line_PRA.
Two findings from doing it, both reported rather than smoothed over:
* io._MODES_ATTR_FRAME_KEYS only knew about "salt_diagnostics", so the new
salt_varying_diagnostics frame fell through to the payload's default=str
fallback and came back from HDF5 as a *string* -- the same class of silent
loss the sidecar was added to fix. Added, with a comment saying a missing
key fails this way rather than loudly.
* Multimode does not converge yet, and the docstring and a runtime warning now
say so. Single-mode is solid (15/15 pumps, residuals < 1e-6); with six
candidates on line_PRA the residual sits near 1.7 and the per-mode
intensities thrash between pumps while their sum still looks monotone --
which is exactly the kind of plausible-looking wrong answer worth refusing
to ship quietly. The cause is the active set, not the operator: a mode is
admitted as soon as D0 passes its *linear* threshold, so modes that should
not lase are handed to the solver and nothing drives them to zero.
_full_salt_newton_impl admits on net gain against the current saturated
background and drops collapsed modes; porting that is the next step.
211 tests pass.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The multimode sweep converged to residuals near 1.7 and per-mode intensities that thrashed between pumps. Two causes, and the second is the interesting one because it reports *success*. 1. The active set was chosen from linear thresholds, handing the solver modes that cannot lase and forcing an impossible lambda_1(k) = 0 on them. Replaced with a self-consistent set: a set is accepted only when every member both lases and satisfies its own residual; failing members are dropped, and candidates are admitted one at a time and kept only if the enlarged set still holds. That alone took the residuals from ~1.7 to ~5e-7. 2. Even then the per-mode curves jumped (one mode reading 0.403, 0.019, 0.392 at consecutive pumps) while the residual stayed at 1e-7 and the summed output stayed monotone. The cause was two active modes drifting onto the *same* k, where the system is degenerate: any split of intensity between them satisfies the equations. The residual cannot see it and neither can the total -- only the per-mode curve shows it. Bounding how far each k may travel fixes it, from the mode spacing, the gain linewidth gamma_perp, and a fraction of k. The same bound is needed for a *single* mode: unbounded, the solve walks to k ~ 0, where the operator is degenerate, and reports a 2.4e-8 residual and converged=True for something that is not a lasing mode at all. Found because switching to bounded least-squares (trf, so a >= 0 is a real bound rather than a clip applied after the solver has already converged) exposed it. line_PRA now gives two lasing modes with strictly monotone L--I curves and residuals of 3e-7..8e-7 at every pump -- the two-mode count being the published Ge-Chong-Stone result. Three regression tests pin the degenerate-k failures, since both of them otherwise look like success. 214 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
line_PRA end to end: 10 m 11.7 s -> 9.1 s. Verified independently of the change:
both L--I curves reproduce the baseline (max rel diff 1.2e-11, bar was 1e-4),
worst residual 8.021e-07 unchanged, 8/8 pumps converged, two lasing modes.
214 tests pass; ruff clean.
Four changes, in descending payoff:
* Vectorised the Magnus propagator and batched it over edges
(edge_transfer_matrices). The sub-interval propagators for all varying edges
become one (n_edges, n_steps, 2, 2) array and the ordered product is n_steps
batched matmuls for the whole graph rather than n_edges * n_steps scalar
ones. The Magnus generator is written in closed form -- for the Helmholtz
system [A1, A2] = (a2 - a1) diag(1, -1), so no matrix commutator is needed.
Bit-identical: same reduction and arithmetic order, relerr 0.0 at n_steps
16/64/256.
* Net-gain admission. A candidate is screened by finding the complex-k root of
lambda_1 on the incumbents' saturated background and admitted only if
alpha < -1e-6, instead of by running a complete trial solve. 35 solves -> 8.
* The per-residual CubicSpline rebuild became a cached resampling matrix:
spline interpolation is linear in the sampled values, and the grid and query
points are fixed while only the values move.
* DENSE_NULL_MAX 400 -> 40, from a measured crossover on the real saturated
operator (n=129: 22.6 ms dense vs 0.98 ms ARPACK), with a seeded start and a
dense fallback for the shift-invert that cannot factorise an exactly singular
operator.
Decomposition, same driver: numerics alone 600 s -> 97 s, active set a further
97 s -> 4.8 s. The 35 -> 8 solve count is only 4.4x but 20x in time, because the
skipped solves were the expensive failing ones -- a non-lasing candidate burns
all 25 outer iterations without converging.
Two real bugs fixed on the way. _smallest_eigenpair's ARPACK branch had no
seeded start vector, so |lambda_1| was not reproducible run to run and the
finite-difference Jacobian inside least_squares read that as noise -- latent
while nothing exceeded 400 nodes, live the moment the threshold dropped. And the
shrink loop ran solve_salt_varying twice with identical arguments, discarding a
solution on failure and re-running the same deterministic solve just to read
argmin(amplitudes); attempt now returns (solution, ok).
Rejected, with numbers, in the code comments: a log-depth pairwise product (2x
at 4 edges, 1.0x at 243 -- the win evaporates exactly at production scale); an
unbounded resampling cache (4 GB and 45 s of thrashing on buffon at n_steps
1024, now budgeted at 64 MB); naive inverse iteration (37x faster, 8% wrong).
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
solve_salt_varying only sees the active set, so with one mode active it has no
spacing to derive a k bound from and falls back to gamma_perp. On a dense
spectrum that is enormous: the buffon's candidates sit 7.6e-4 apart against a
0.5 fallback, 658x the spacing, so a single active mode wanders across every
candidate.
A single global cap from the *minimum* spacing is the obvious fix and is also
wrong -- it punishes well-separated modes for a near-degenerate pair elsewhere.
On the buffon two candidates are 7.6e-4 apart while the others are ~20x further
out, and a global 0.2*min-spacing cap (1.5e-4) left nothing able to lase at all.
So k_window_cap now accepts one value per mode, and the sweep derives each from
that mode's own nearest neighbour among all candidates.
line_PRA is unaffected: same L--I curves, same 8.021e-07 worst residual, 8/8
converged. 214 tests pass.
This does NOT fix the buffon, and the honest diagnosis is now sharper. With
correct per-mode bounds and a 20-point pump grid the solver still lands on
spurious large-amplitude roots from the first pump onward:
D0 = 1.02x threshold a = 6.29 res = 8.8e-07 converged
D0 = 1.26x a = 137.79 res = 4.7e-07 converged
D0 = 2.09x a = 418.26 res = 2.8e-01 not converged
At 1.02x threshold the physical amplitude must approach zero; 6.29 is a
different root of the same equations, and it satisfies them to 1e-7. Refining
the pump grid makes it worse, not better, so this is not continuation
coarseness -- the initial amplitude guess (a flat 1e-2) sits in the wrong basin
on this graph, where the physical near-threshold amplitude is much smaller and a
spurious branch is nearby. The next step is to seed `a` from the linear
prediction (D0/D0_thr - 1)/T_mumu as modes._full_salt_newton_impl's
onset_amplitude does, rather than from a constant.
Also noted while reproducing: graph.graph["params"]["pump"] is left unmasked on
boundary edges (26 of the buffon's 243), so a caller reading it directly instead
of the pipeline's pump profile hands a varying profile to a lead and trips
construct_laplacian_varying's boundary guard. The pipeline itself is correct.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
… one branch Two defects found by auditing the operator at 208 nodes, where everything to date had only been validated at 8. 1. x_scale. The residual's sensitivity to k and to a differ by ~6e5 on the buffon (|dlam/dk| = 8.2e2 against |dlam/da| = 1.3e-3), because the amplitude's unit is set by the pump-region norm: mean |E|^2 is 5.4e-4 over 2500 units of pumped length there against 1.04 over 1.0 on line_PRA, so the natural amplitude is ~1900x larger. least_squares ran with the default isotropic x_scale=1.0, taking steps sized for k that are useless for a, and the amplitude never left its initial guess. This inverts the earlier reading of the failure: a = 6.29 at 1.02x threshold was not a physically impossible large amplitude, it was an unmoved small one -- 6.29 is 0.0033 in line_PRA units. Amplitudes on this graph are *supposed* to be in the hundreds, and any intuition calibrated on line_PRA will misread them. 2. Branch bound. _smallest_eigenpair returns argmin |lambda|, a min over analytic branches, so past the k where two branches cross, the eigenvalue being root-found belongs to the neighbouring mode and is discontinuous there. The crossing was measured at 0.395 of the candidate spacing, so the 0.4 bound put it 3.1e-5 *inside* the box: the solve slid onto the neighbour (eigenvector overlap 0.89 with the wrong mode), pinned at the bound and drove the amplitude to zero. Tightened to 0.25. line_PRA is unchanged by both: same L--I curves, same 8.021e-07 worst residual, 8/8 converged, two lasing modes. 214 tests pass. x_scale="jac" costs Jacobian evaluations, 9.1 s -> 16.9 s on that example. Ruled out by the same audit, and worth not re-investigating: the ARPACK shift-invert path agrees with dense eig to 6.6e-14 in eigenvalue and 2.2e-16 in eigenvector over 278 calls of a real buffon solve, with the dense fallback never taken; the M12 Dirichlet pole is *better* conditioned on the buffon than on line_PRA (min |M12| 6.5e-4 against 4.2e-4, and the smallest are unpumped boundary edges); and the unconjugated pump norm does not suffer cancellation there (int|psi|^2 / |int psi^2| = 1.24-1.29, same order as line_PRA). Not yet fixed. The buffon still does not produce a usable L--I curve: at 1.02x threshold it now converges to a = 1.69 at residual 4.6e-7, while an independent 2-variable root find on the same frozen field reaches a = 343 -- two distinct roots, and which is the physical branch is unresolved. n_steps=128 is also under- resolved there (2.3e-3 relative in lambda_1 against 3.2e-9 on line_PRA, ~1.1% in the amplitude), so 256 is the honest minimum for that graph. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…velength
Ground truth for the buffon's lowest mode was established independently
(amplitude continuation from a = 0, cross-checked against a root map and against
the analytic onset slope 1/T_00; roots held to |lambda_1| ~ 1e-13, three orders
tighter than the module's target). The correct curve is
D0/D0_thr 1.02 1.05 1.10 1.26 2.09
a 1.431 3.092 5.623 15.558 68.077
monotone, bending below the linear model, meeting its slope at threshold to
0.8%. Every amplitude this solver has reported was wrong -- and so was my
reading of them: 6.29 at 1.02x was not "physically impossible", it was 4.4x too
large, and the earlier claim that a ~ 343 was the root is also wrong.
Two independently measured defects fixed:
* The k bound came from the *candidate list*, which is only as complete as the
mode search behind it. On this graph it is not remotely complete: Weyl gives
~120 modes in k = [10.63, 10.73] (mean spacing pi/L_opt = 8.4e-4) where the
fixture carries 4, so a bound keyed off those 4 still admitted whole
unlisted modes and the solver converged to genuine 1e-7 roots belonging to
other modes. Now min(Weyl mean spacing, nearest listed candidate). A
residual pins k to 2.5e-10 here and says nothing about being on the intended
branch, which is why this was invisible.
* n_steps is a per-edge sample count, not a resolution, and the graphs differ
by two orders of magnitude in wavelengths per edge: 0.92 on line_PRA against
50.9 on the buffon, so the same default is 69 samples per wavelength on one
and 1.3 on the other. Measured against the onset slope, that is a 47%
amplitude error at n_steps=64, 0.2% at 256. _resolved_n_steps now raises
n_steps until a pumped edge gets 8 samples per wavelength, bounded by a cap.
line_PRA needs 8, far below its default, so it is untouched.
Together these take the buffon from wrong at every pump to correct within 1% at
the two pumps nearest threshold:
1.02x a = 1.420 (truth 1.431, -0.8%) residual 4.6e-07 converged
1.05x a = 3.065 (truth 3.092, -0.9%) residual 6.5e-07 converged
Still broken from 1.10x up, and the cause is now narrowed: the amplitude jumps
to 52.7 where truth is 5.62, which is neither resolution nor the k window (both
are fixed above and the failure is unchanged by either). It is a branch jump in
a. The next suspect is the seed field -- solve_salt_varying seeds the hole
burning from node_solution_varying(k, graph, None), the *unpumped* operator,
where |lambda_1| = 4.57 and the state is not a mode at all; its T_eff is 2.6x
too small.
line_PRA unchanged throughout: same L--I curves, same 8.021e-07 worst residual,
8/8 converged. 214 tests pass.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The buffon's L-I curve above 1.05x threshold is wrong: the solve reports
a = 55.3 where an independent amplitude continuation gives 5.62. This
lands one measured part of the cause, and records what the same
measurements rule out.
What the measurements say, on the production buffon at D0 = 1.10x
threshold, started from the converged 1.05x state (a = 3.06):
* The equations are right. Fixing `a` and solving for (k, D0) --
the same two equations Re/Im lambda_1 = 0, the same `_lam_varying`
residual, only a different pair of unknowns held -- reproduces the
whole L-I curve to ~2%: a = 1.436, 3.22, 5.5, 15.47, 66.8 against a
truth of 1.431, 3.092, 5.623, 15.558, 68.077 at 1.02, 1.05, 1.10,
1.26, 2.09x threshold, with a near-threshold slope of 75.5 against
the analytic 1/T00 = 74.01.
* The frozen field is not at fault. Holding k, the frozen-field
residual has a clean minimum exactly at the true a = 5.62
(|lambda| 0.287 -> 0.018 -> 0.356 at a = 3.06, 5.62, 10), so the
root is present and correctly placed in the problem the solver is
actually given.
* Past a ~ 12 the residual flattens at |lambda| ~ 0.47 while
|dlambda/da| collapses from 9.6e-2 to 3.4e-3, a 28x drop.
Hence the ceiling: an amplitude may rise by at most 4x per frozen-field
solve, expanded only while expanding actually lowers the least-squares
cost. Expanding merely because the solution sits on the bound is what
the runaway wants -- the flat region is above every finite ceiling, so
"pinned" stays true all the way up.
This is necessary but not sufficient, and the same runs say why: with
the ceiling in place the buffon still reaches a = 52.7, coming back with
k pinned at *its* cap (+2.57e-4 against a true +7.3e-5). The runaway is
a joint (k, a) direction, not an amplitude excursion; the residual scan
sees a clean minimum precisely because it holds k fixed. Closing that
needs the solve to start inside the right basin, which only the caller
can arrange.
line_PRA is unchanged to every printed digit (0.16499..0.90914 and
0.05459..0.92138, 8/8 pumps converged, worst residual 8.02e-07, same
iteration counts 13/14/15/15/17/17/17); it never expands the ceiling.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…seed
Two changes, both aimed at the buffon L-I curve above 1.05x threshold.
It is better and still wrong; the numbers and the new open question are
below.
1. `_predicted_amplitude` seeds each pump by extrapolating the (D0, a)
pairs the sweep has already accepted, instead of reusing the previous
pump's amplitude. From the 1.02x and 1.05x solutions (a = 1.419,
3.063) the prediction for 1.10x is 5.80 against a true 5.623.
2. A frozen-field solve that returns `k` sitting on its own cap is no
longer accepted, whatever its cost. The cap marks the edge of the
mode's analytic branch, so a k pinned there means least_squares
wanted to leave the branch and was merely stopped -- it reports the
boundary, not a root. When every solve in a round is such a corner,
(k, a) are held and only the field is refreshed, since a stale field
is precisely why no interior root was available.
Effect at D0 = 1.10x threshold on the production buffon:
before a = 52.7 k-k0 = +2.57e-4 (on the cap) res 4.3e-01 no
after a = 22.7 k-k0 = +4.45e-5 (interior) res 6.8e-07 yes
So the runaway and the pinning are gone. What replaced them is worse to
diagnose and better to have: a *converged* answer, at 6.8e-07, that is
still not the physical one (truth 5.623).
Note (2) alone did this; (1) did not move the result. Seeding the solve
at a = 5.80 -- inside the right basin, from the converged 1.05x field --
still left for a = 52.7 before this change. That rules out the initial
guess as the cause and is why the prediction is kept only as the cheap
improvement it is, not as the fix.
Open question this exposes, for the next commit: whether a = 22.7 is a
genuine second root of the self-consistent equations at this pump, or an
artifact of the amplitude ceiling -- 22.732 against a ceiling of 4 x
5.801 = 23.204 is 98% of it, which is close enough to demand a check.
If it is a real second root, no amount of local-solver hygiene picks the
right one and the physical branch has to be selected by continuity from
a = 0, which is what the independent amplitude continuation does.
line_PRA is unchanged to every printed digit (0.16499..0.90914 and
0.05459..0.92138, 8/8 converged, worst residual 8.02e-07, iteration
counts 13/14/15/15/17/17/17), 6.1 s. 210 unit tests pass.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Records what the buffon failure is, with a reproducer, after three
plausible explanations were measured and ruled out.
The finding: fix `a` and solve for (k, D0) instead of fixing D0 and
solving for (k, a) -- same two equations, same `_lam_varying` residual,
different pair of unknowns held -- and the resulting D0(a) is STRICTLY
INCREASING at all 56 continuation points from a = 0.2 to a = 70. A
strictly monotone D0(a) is a bijection, so each pump has exactly one
amplitude on this branch. The continuation places the solver's converged
a = 22.73 at D0 = 1.38x threshold, not the 1.10x it was asked for.
So the solver returns the amplitude belonging to a different pump, and
its residual of 6.8e-07 says nothing: the 1.10x row is converged and
wrong by 300%. That makes this a solver problem rather than an equations
problem, and it is why the earlier "converged" reading was not evidence.
The continuation is independently anchored: near-threshold slope 75.5
against the analytic 1/T00 = 74.01, and max |k - k0| = 7.3e-05 over the
whole trace, inside the 2.09e-04 k cap, so it never leaves the mode.
Also recorded, because both produced confident wrong numbers here and
cost several probe rewrites:
* a probe that lets k travel silently measures the NEIGHBOURING mode --
at a = 20 there is a genuine root at k - k0 = +6.5e-04, three times
the cap, and two probe versions tracked it and reported smooth,
monotone, meaningless curves;
* scanning |lambda| at real k cannot answer a branch question, because
for an `a` that is not a solution there is no real-k root and the
number reported is the bound k ran into.
Adds:
* examples/buffon/buffon_narrow -- the narrow-k-window buffon variant
the measurements are made on, as a config rather than a description;
* examples/audit/probe_varying_amplitude_branch.py -- traces the
branch, checks D0(a) is monotone, and diffs the solver against it;
* AUDIT.md section 10 and the audit README entry.
No library change. line_PRA remains correct and unchanged throughout.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The fixed-D0 solve selects the wrong branch above ~1.05x threshold. It is
not an equations problem -- AUDIT.md section 10 shows D0(a) is strictly
monotone along the physical branch, so each pump has exactly one amplitude
on it, and the solver was returning the amplitude belonging to a different
pump (a = 22.73 at 1.10x, which the branch places at 1.38x). Selecting the
right one requires continuity from a = 0, which a fixed-D0 solve cannot
express.
So invert which unknowns are held. The unknown vector becomes
x = (k_0 ... k_{M-1}, u_1 ... u_{M-1}, D0)
-- still 2M unknowns for 2M equations -- with the amplitudes recovered from
a gauge-fixed weight vector, a_mu = s * u_mu / sum(u) and u_g = 1, so
sum(a) = s by construction: no penalty term and no constraint equation to
weight against the residuals. `s` is the continuation parameter, D0 is
solved *for*, and a secant on `s` walks the achieved D0 onto the requested
one.
Only the TOTAL scale is gauge-fixed, never one mode's amplitude. Pinning a
single mode is the seemingly equivalent thing and it is not: it left
line_PRA's first five pumps byte-identical and then collapsed mode 1 to
zero and invented a third mode. Every non-gauge mode can still reach a = 0
and be rejected by the active-set logic, since u_mu = 0 is inside bounds.
Convergence stays judged by salt_residuals_varying at the REQUESTED D0, so
a solve whose secant has not landed reports a large residual and
converged = False rather than a small residual at a pump nobody asked for.
`outer` 25 -> 40: the continuation needs more field refreshes than the
fixed-D0 solve, and the buffon's 1.05x pump was consuming exactly 25.
Measured on the buffon (examples/buffon/buffon_narrow, n_steps = 512)
against the independent amplitude continuation, all converged:
D0/D0_thr a truth err k - k0 residual
1.02 1.419 1.431 -0.8% +1.59e-05 9.2e-08
1.05 3.063 3.092 -1.0% +4.76e-05 2.8e-07
1.10 5.576 5.623 -0.8% +7.32e-05 4.7e-07
1.26 15.388 15.558 -1.1% +2.82e-05 8.2e-07
1.10x was previously 22.73 (+304%) and 1.26x was 85.67 (+451%, not
converged). Every k - k0 stays well inside the 2.09e-04 cap, and the 1.10x
value of +7.32e-05 reproduces the independently derived +7.3e-05 -- the
check that matters, since the earlier wrong answers all came back with k
pinned on the cap.
The 2.09x pump (truth 68.077, a 4.4x amplitude step) had not finished at
commit time and is not yet verified.
line_PRA is unchanged to every printed digit -- two lasing modes,
0.16499..0.90914 and 0.05459..0.92138, 8/8 pumps converged, worst residual
9.16e-07 -- at 9.8 s against 6.4 s. 210 unit tests pass. Buffon cost is
~2 min per pump against ~1 min.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
…dependent solver
Two things, both prompted by a fair challenge to how "ground truth" was
being used.
CORRECTION. The buffon reference curve (a = 1.431, 3.092, 5.623, 15.558,
68.077) was produced by amplitude continuation through netsalt's OWN
varying operator -- same construct_laplacian_varying, same
saturated_eps_profiles, same _lam_varying. Calling it ground truth
overstated it. It independently checks which solution BRANCH is selected,
which was exactly the defect, so a 300% disagreement between two
parametrisations of one model was a genuine contradiction. It does not
check the model: a systematic operator error would be shared and the two
would agree perfectly. The docstring and AUDIT.md section 10 now say so,
and name the real anchors -- the linear competition matrix's onset slope
near threshold (1/T00 = 74.01 vs 75.5, a different code path) and
independent_salt/ for the model itself.
Also corrected: both tables still showed the PRE-fix solver numbers, and
neither recorded that the 2.09x pump still fails. It is a 4.4x amplitude
step, too large for the secant to bridge; the solve holds at the incoming
value and reports residual 1.2 with converged = False -- honest failure,
not a wrong answer -- after exhausting 40 outer iterations in 38 minutes.
The working range is ~1.26x threshold, against ~1.05x before, not 2x.
NEW TOOLING. examples/audit/independent_salt/ validates netsalt against a
solver sharing no code with it (numpy and scipy only), but only ever the
OVERSAMPLED path -- absolute intensities to 1.6e-4 above threshold over
122 pump points. The varying path has never been run against it, which
AUDIT.md section 9 flags as outstanding. Adds:
* step5v_netsalt_varying_run.py -- the varying-operator sweep on the
same Fabry-Perot, emitting step 5's schema. The physical intensity is
a * int_cavity f(x) dx, by Simpson quadrature over each edge's sample
grid rather than l_e * f_e, since here f is sampled ALONG the edge and
an O(h^2) rule would discard the fourth-order propagator's accuracy;
* step8v_compare_varying.py -- direct comparison, no Richardson stage,
because the varying path's convergence parameter is n_steps rather
than an oversampling resolution.
Generated results_*.json and *.h5 under independent_salt/ are gitignored,
matching how the other steps' artefacts are treated.
No library change.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Section 10's reference curve ran through netsalt's own varying operator, so
it checked which solution branch was selected and not the model itself. This
runs the varying path against examples/audit/independent_salt/ -- a solver
written from the equations, sharing no code with netsalt, verified to 0-4e-15
against the closed-form Fabry-Perot spectrum -- which had only ever been used
on the oversampled path.
Fabry-Perot, n_steps = 64, D0 = 0.58 .. 1.40, two co-lasing modes, 72
(pump, mode) comparisons:
lasing frequency k_mu : median 4.4e-06 max 5.7e-06
modal intensity I_mu : median 1.1e-04 max 7.8e-03
Absolute intensities, not ratios -- int|Psi_mu|^2 dx on both sides, so SALT's
denominator fixes the scale and no normalisation is free. At D0 = 0.70 the
varying path gives 7.804615e-2 against the independent solver's 7.804230e-2.
The oversampled path's published figures on the same case are median 1.6e-4 /
max 8.1e-3, so the varying operator matches independent truth at least as
well on a matrix that never grows. The worst point (D0 = 0.82, mode 2) sits
in the pump steps straddling the second mode's turn-on where its intensity is
near zero -- where the oversampled path's outliers also sit.
One bug in the comparison, worth naming because it failed silently: the
independent solver normalises int|phi|^2 = 1, so its `a_mu` IS the physical
intensity int|Psi_mu|^2 dx, while this side reports that as `I_mu` and keeps
`a_raw` for the un-normalised amplitude. Looking up `I_mu` on the reference
matched nothing and the script reported "no overlapping pump points" rather
than a wrong number -- but comparing `a` to `a` would have compared two
different normalisation conventions and produced a plausible, wrong answer.
This validates the varying MODEL above threshold and multimode on a 1D
Fabry-Perot. It does not independently validate the buffon, whose geometry is
far harsher. What it does mean is that section 10's continuation reference no
longer rests only on self-consistency.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The continuation walks the amplitude and reads off the pump it lases at, so
the distance it must travel is set by how far D0_target sits from the pump
the caller's seed actually solves. On the buffon a 1.26x -> 2.09x request is
a 4.4x climb in amplitude (15.4 -> 68) from a seed whose own pump is 0.60 of
the target, and the secant does not bridge it in 40 outer iterations: it
holds at the seed and reports converged = False.
Splitting the interval geometrically -- geometric because the amplitude
climbs multiplicatively, so a 4.4x step halves into two of 2.1x -- turns one
step the secant cannot take into two ordinary ones. The retry runs only when
the solve failed AND the pump ratio exceeds _SUBSTEP_MIN_RATIO, recurses at
most _SUBSTEP_MAX_DEPTH times (16 sub-solves), and keeps the original answer
if the sub-stepped one is no better, so it can only help.
The first bracket is reused as the seed's own pump: _bracket_D0 already
answers "at what pump does this amplitude lase?" on the first outer
iteration, which is exactly the interval's lower end.
Verified NOT to regress:
* line_PRA byte-identical -- 0.16499 .. 0.90914 and 0.05459 .. 0.92138,
8/8 pumps converged, 7.2 s;
* 210 unit tests pass;
* buffon 1.02x/1.05x/1.10x/1.26x unchanged at -0.8%/-1.0%/-0.8%/-1.1%
against the continuation branch.
NOT yet verified: the 2.09x pump this exists for. That solve is still
running at the time of this commit and is committed only so the work is not
lost to a container reclaim -- the sweep has been killed three times. If it
does not hold, this is a revert, not a fix.
Refs #52, #53.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
It does not fix what it was written for, and it costs time to fail. On the buffon's 2.09x pump the retry fired, the mid-point solve failed too, and the recursion took 83 minutes to reach the same answer the single attempt reached in 38: a = 15.388 held at the seed, residual 1.2, converged = False. So the 2.09x failure is NOT a step-length problem, which is what sub-stepping assumes. The evidence against that reading was already available and I did not weigh it: the 1.10x -> 1.26x step is a 2.8x amplitude climb and succeeds, while the mid-point this retry chose (~1.62x, a 2.3x climb from the same seed) fails. A smaller step failing where a larger one succeeded rules out step length. What the continuation table suggests instead: between a = 15 and a = 68 the frequency pull turns around -- k - k0 runs +3.2e-05 at a = 15, down to -7.0e-05 near a = 44, back to -2.6e-05 at a = 68. The secant assumes D0(s) is smooth and monotone, which it is, but the *field* has to follow a k that reverses direction, and the frozen-field refresh may not track that. That is the next thing to test, and it is a different mechanism. Reverting rather than leaving it disabled: it is dead weight on every failing solve, and the diagnosis it encodes is wrong. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
This is what the buffon's large-pump failure actually was, and it was never a step-length problem -- the sub-stepping reverted in 9893bdd assumed it was. The field refresh is a fixed-point iteration on the hole-burning field, and at high amplitude the shipped damping = 0.7 makes it DIVERGE. Instrumented on the buffon at 2.09x threshold, seeded from the converged 1.26x state, the relative field movement pins at 6.618e-01 and holds there to four digits for every one of the 40 outer iterations, against a 1e-3 tolerance: the field is replaced wholesale each round instead of settling. Since the scale may only step once the field has caught up, the continuation stayed frozen at its seed -- the debug trace shows s = 15.402 unchanged while the inner solve converged beautifully (cost 1e-26, off its k bound) to the pump that amplitude really lases at, 0.51 D0_target. Nothing was wrong with the solve; the gate in front of it never opened. Under-relaxing fixes it outright. Halving on demand rather than lowering the default keeps the cases where 0.7 already converges untouched, and easing it back up after a refresh that settles at once stops the backoff ratcheting -- without recovery, one hard state permanently slowed every later one and left 1.26x (21 iterations at the shipped damping) burning all 40 and reporting failure with the right answer. `outer` 40 -> 80: with under-relaxation the continuation genuinely needs more rounds, and 2.09x now uses 35. Buffon, n_steps = 512, against the continuation branch -- every pump converged, where 1.10x upward previously did not: D0/thr a truth err residual 1.02 1.419 1.431 -0.8% 9.2e-08 1.05 3.063 3.092 -1.0% 2.8e-07 1.10 5.576 5.623 -0.8% 4.7e-07 1.26 15.387 15.558 -1.1% 1.0e-07 2.09 66.781 68.077 -1.9% 4.1e-07 Worst error 1.9%, down from 77.4%. The 2.09x row also lands k - k0 = -3.2e-05 against the branch's -2.6e-05, i.e. on the far side of the frequency-pull reversal (the pull peaks at +7.3e-05 near a = 5.6 and comes back down), so it is following the branch and not merely reaching a similar amplitude. line_PRA byte-identical: 0.16499 .. 0.90914 and 0.05459 .. 0.92138, 8/8 pumps converged, 13-22 iterations, 8.1 s. 210 unit tests pass. Refs #52, #53. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
The docstring still said "still fails" and carried the pre-fix numbers. Updates the table (every pump now converges, worst error 1.9%, residuals 9.2e-08 .. 4.7e-07) and replaces the step-length explanation, which was wrong. Sub-stepping through the geometric midpoint failed while the 2.8x step from 1.10x to 1.26x had always succeeded -- a smaller step failing where a larger one succeeds rules out step length. The cause was the field refresh diverging at high amplitude: relative movement pinned at 6.6e-01 against a 1e-3 tolerance, so the gate that lets the scale step never opened and the continuation stayed frozen at its seed, while the inner solve converged perfectly to the pump that seed amplitude really lases at. Also records why the solved k - k0 = -3.21e-05 at 2.09x is evidence rather than coincidence: the frequency pull reverses along this branch, peaking at +7.3e-05 near a = 5.6 before coming back down, so landing on the far side of that turn means the solver followed the branch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
`relax` is what distinguishes "the field refresh is backing off" from "the field refresh is fine and the cost is elsewhere", and without it in the trace that question needs a code edit to answer. It settled one hypothesis already: multimode solves got slower after the under-relaxation landed, and the obvious suspect was the backoff firing where it was not needed -- M = 1..4 all converged at that pump before under-relaxation existed. The trace shows relax = 0.7000 throughout with dfield = 5e-04, comfortably inside the gate, so the backoff never fires there and the multimode cost is intrinsic to the solve rather than a regression from that fix. Debug-branch only; no behaviour change. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Speed-up audit. Profiling a 3-mode solve put 83% of the time in building the
operator and ~1% in the eigensolve, with 69% cumulative inside
`_SaturatedEdgeProfile.__call__` -- the cubic interpolation of D0_eff from the
uniform field-sample grid onto the Magnus abscissae. So the eigenproblem is not
the cost and an analytic Jacobian would have optimised the 1%.
The access pattern is what makes this cacheable. `saturated_eps_profiles`
builds these objects once per residual evaluation; `_lam_varying` then rebinds
`gain` and rebuilds the operator once per mode. Each object is therefore called
2M times per residual (magnus4 asks for two Gauss point sets) against only TWO
distinct query arrays, alternating between them. `_d0_eff` is fixed for the
object's life, so all M repeats interpolate identical values.
Caching per object collapses that to two interpolations. Measured on the
production-window buffon (208 nodes, 217 pumped edges, n_steps resolved to
408), M = 3, outer = 4:
CubicSpline constructions 273470 -> 83750 (3.3x, = M as predicted)
wall 109.2s -> 93.6s (1.17x)
Amplitudes are bit-identical (0.012538578 and the two zero modes to 8 digits),
line_PRA is unchanged to every digit, 210 unit tests pass.
Note a single-entry cache is worth exactly nothing here, and an earlier attempt
measured precisely 1.00x: the two Gauss sets alternate, so a one-slot memo
misses on every call. The dict is the point.
Also measured and NOT taken:
* raising CUBIC_RESAMPLE_CACHE_BYTES 64 MiB -> 1 GiB. At the resolved
n_steps = 408 an entry is ~1.3 MiB and ~434 are wanted, so only ~48 fit and
the rest miss forever. Fixing that is worth a further 1.075x for 2.5x the
peak RSS (0.25 -> 0.63 GiB), which is not a trade worth making by default.
* replacing the dense resampling matmul with a local B-spline evaluated per
call: 0.56x, i.e. slower. The cached matmul avoids rebuilding the spline,
which is exactly what the local form has to do.
* SALT_VARYING_SAMPLES_PER_WAVELENGTH 8 -> 5: 1.31x, and the module's own
convergence table puts the cost under 1% (slope 74.19 against the analytic
74.01). Left alone because it trades accuracy for speed and is already
exposed as `intensity_varying_samples_per_wavelength`.
The ceiling, measured by handing the propagator values already on its abscissae:
the operator build is 12.5x faster (0.659s -> 0.053s), ~2.74x overall. Reaching
it means sampling the field where the propagator asks rather than on a uniform
grid and interpolating -- which also moves the pump-region norm from Simpson to
Gauss quadrature. That is an architectural change, recorded here rather than
attempted.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Takes the ceiling the speed-up audit measured. 69% of a multimode solve was
`_SaturatedEdgeProfile.__call__`, cubic-interpolating the saturated D0_eff from
the uniform field-sample grid onto the Magnus abscissae -- a mismatch that
exists only because `edge_field_samples` returns psi at the sub-interval
boundaries while the propagator evaluates eps at Gauss points *inside* each
sub-interval.
So sample |E|^2 at the abscissae in the first place. The interpolation moves
from the profile callable to `edge_field_profiles`, i.e. from once per edge per
mode per residual evaluation to once per edge per field refresh, and the
callable becomes a slice lookup. Handing the propagator values already on its
abscissae was measured at 12.5x on the operator build alone (0.659s -> 0.053s).
Measured end to end on the production-window buffon (208 nodes, 217 pumped
edges, n_steps resolved to 408), M = 3, outer = 4:
109.2s -> 44.0s (2.48x; 93.6s with the per-object cache alone)
and on the buffon truth ladder at n_steps = 512, 878s -> 461s (1.9x) with the
answers unchanged: -0.8, -1.0, -0.8, -1.1, -1.9 % against the continuation
branch, every pump converged, worst residual 7.7e-07.
This is not bit-identical, and it should not be. The interpolation is the same
order but no longer the same operation: |E|^2 is resampled and the saturation
formed pointwise at the abscissae, where before the *result* of the saturation
was resampled. line_PRA's intensities are unchanged to every printed digit
(0.16499..0.90914 and 0.05459..0.92138, same iteration counts, 8/8 converged)
and its residuals move in the 7th digit. The buffon's lasing amplitude holds to
8 digits; its two non-lasing modes move from ~1e-21 to ~1e-17, both numerical
zero.
A fallback interpolation is kept for callers whose profiles came from a
different n_steps or method, so mixing resolutions stays correct rather than
merely fast.
Two unit tests needed their quadrature corrected, not their expectation: both
summarised a profile with a uniform-grid rule (`trapezoid`/`simpson` at
dx = L/(n-1)) that is simply the wrong rule on Gauss abscissae. Two-point Gauss
weights every sample equally, so the edge mean is the arithmetic mean. Checked
before touching them -- the normalisation invariant they exist to protect still
holds to 2.2e-08 against `modes._single_mode_field_intensity`.
210 unit tests pass.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
Issue #54 asks whether full SALT handles three or more co-lasing modes, and nothing in the repo had reached it: line_PRA lases two by construction, and the independent Fabry-Perot cannot reach three at all (uniform index and pump make the competition nearly rank-1). The buffon's spectrum is dense enough to ask -- all twelve candidate thresholds fall within 4.7% of each other, so competition rather than threshold ordering decides which modes lase. Six modes lase simultaneously at 1.05x the lowest threshold, converged to 5.0e-07 in 17 iterations, and agree with the linear competition matrix -- a different code path, no transfer matrices, no varying operator -- to within +-7%. SALT sits above linear on the strong modes and below on the weak ones, which is the sign the competition correction should have. Also records that the odd-looking intensity ordering (mode 1 carrying ~9x mode 0 while sitting closer to its own threshold) appears in BOTH solvers, so it is pump overlap rather than a solver artefact; and that seeding the set from the linear model is what makes this affordable -- the same six modes asked for at 1.01x, where only one of them lases, had not finished in 14 minutes against 6 for the well-posed pump. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015cB7hwiL9WmKk4nH6rULvS
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Groundwork for #52 / #53. Small and self-contained: one new module, its tests, and one audit reproducer. Nothing is wired into the solver yet — that is the follow-up, and this PR is what establishes the design for it.
A correction first
My earlier comment on #52 recommended computing the within-edge hole-burning integral analytically per edge, "as
compute_mode_competition_matrixalready does". That was wrong, and it would have sent the follow-up work down a dead end. There is no such integral in the operator:One
k_eper edge. That closed form is the exact solution of the 1D Helmholtz equation only when ε is constant along the edge — it is precisely what makes the quantum-graph secular matrix exact.dispersion_relation_pump_saturatedlikewise takesD0_effas one scalar per edge.So the real situation is structural, not a badly-chosen budget: spatial hole burning makes ε vary within an edge, breaking the piecewise-constant assumption the whole method rests on.
oversample_graphrestores piecewise-constancy by subdividing, and pays for it in the size of the eigenproblem — 76803 nodes against 243 edges on the production buffon, which is what puts full SALT out of reach there.(The competition matrix's analytic integral is a different object: a known field product over an edge at fixed ε. That factorises; this does not.)
What this adds
netsalt/edge_propagator.py— the 2×2 transfer matrix of an edge whose permittivity varies along it. An edge with varying ε still has an exact propagator; it is simply not available in closed form. Computing it per edge keeps the eigenproblem at its original size and makes the sub-interval count local, independent and parallel.propagator_constant_eps(q, length)— the closed form, in transfer-matrix form.edge_transfer_matrix(k, length, eps, n_steps, method)—epsmay be a scalar or a callable. A scalar is propagated exactly in one step, no discretisation.Two measured results
examples/audit/probe_edge_propagator.py, against a DOP853 reference at rtol 1e-13, on a buffon-like edge (length 11,k = 10.7,n = 1.5→ ~28 oscillations, 4 % saturation ripple):exp(h·A(x_mid))for this system — they agree to round-off, which is what makes the comparison fair. It also means the current scheme is O(h²) and no tuning of the node budget can improve that order.The larger point is where the count lands. Today it multiplies the graph: 243 buffon edges × N sub-edges is the eigenproblem size, and N ≈ 316 already means 76803 nodes. As a per-edge transfer matrix the eigenproblem stays 243 edges.
Tests
Six new tests, on the properties that actually matter:
epsbypasses discretisation entirely;det = 1for any lossless profile;190 tests pass;
ruff checkandruff format --checkclean.Follow-up
Wiring this into
construct_incidence_matrixandconstruct_weight_matrix— where theopen,directedanddirected_reversedboundary models each need handling, and where the passive path must stay byte-identical. Kept separate deliberately: that change carries regression risk on the passive solver, and this one carries none.Generated by Claude Code