diff --git a/.gitignore b/.gitignore index cfe3ff29..6c1dd0b5 100644 --- a/.gitignore +++ b/.gitignore @@ -114,3 +114,12 @@ venv.bak/ *.swo *.png + +# uv lockfile (not pinned for this project) +uv.lock + +# Example run artifacts (generated by create_graph.py / the pipeline) +examples/**/out/ +examples/line_PRA/graph.json +examples/line_PRA/index.yaml +examples/line_PRA/pump.yaml diff --git a/CLAUDE.md b/CLAUDE.md index be9462e4..30e1785b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -20,7 +20,18 @@ mapping. - `modes.py` — mode search driver (`scan_frequencies`, `find_modes`, `find_threshold_lasing_modes`, `pump_trajectories`, `compute_mode_competition_matrix`, `compute_modal_intensities`). Runs - `multiprocessing.Pool` over the scan grid. + `multiprocessing.Pool` over the scan grid. The lasing L–I curves can be + computed by two `intensity_method`s (issue #42, see `doc/source/lasing.rst`): + `linear` (default, the fast near-threshold competition-matrix model) and + `full_salt_newton` + (operator-level nonlinear SALT — solves `(k_μ,a_μ)` so `L_sat` is singular at + real `k_μ`; reduces to `linear` near threshold and bends the curves above it. + Auto-oversamples the within-edge hole burning — the per-edge mean over-clamps; + validated against Ge-Chong-Stone PRA 82, 063824 Eq. 28 on `line_PRA`). The + oversampling is bounded by `oversample_node_cap` (default 3000), so newton + **also runs on the real buffon** (~12 s capped, ~6 min uncapped at λ/4 ≈30k + nodes, lasing its co-lasing modes) — raise the cap / `oversample_resolution` + for accuracy at higher cost. `linear` is still the cheap first pass for counts. - `algorithm.py` — rough mode detection (skimage `peak_local_max`) and two refinement algorithms: `refine_mode_root` (MINPACK ``hybr``, default) and `refine_mode_brownian_ratchet` (legacy random-walk @@ -60,8 +71,11 @@ mapping. `compute_lasing_modes` on a line graph and diffs `out/` against `tests/data/run_simple/out/` with `dir_content_diff` - `examples/` — ready-to-run YAML configs (buffon, ring, wheel, directed, - line_PRA, transfer). Buffon variants use `defaults:` to inherit shared - base configs. + line_PRA, transfer; buffon variants use `defaults:` to inherit shared base + configs) plus script-based per-graph solver comparisons (line_fabry_perot, + ring_leads, tree, two_ring, chaotic_ring, dense_ring; shared helpers in + `examples/_common.py`, figures gitignored and reproduced by each folder's + `run.sh`). - `doc/` — Sphinx sources, published at https://arnaudon.github.io/netSALT/ ## Running things @@ -209,6 +223,27 @@ they are what an "old code" most needs before further work lands on top. - Compute-core tests are still fairly shallow. Adding a test for ``compute_mode_competition_matrix`` and ``find_threshold_lasing_modes`` on a tiny analytic graph would give more regression coverage. +- ~~**Full SALT beyond the linearised competition matrix (issue #42).**~~ + **Landed.** The operator-level ``full_salt_newton`` solver was added next to + the original ``linear`` model: it + solves the real nonlinear SALT eigenproblem (saturated dispersion + ``dispersion_relation_pump_saturated`` + a frozen-field trust-region ``(k,a)`` + solve with a self-consistent active set). All reduce to ``linear`` near + threshold. ``full_salt_newton`` is **validated against Ge-Chong-Stone PRA 82, + 063824 (Eq. 28)** on ``line_PRA``: both lase two modes, intensities matching to a + few % near threshold, with the full-SALT bend-over above it. Crucial gotcha: the + operator-level hole burning must **resolve the within-edge field** — the per-edge + mean over-clamps and spuriously suppresses co-lasing modes (it lased one mode on + ``line_PRA`` until ``oversample_size`` auto-defaulted to a wavelength-resolving + size, ``_auto_oversample_size``). ``benchmark/bench_salt.py`` compares the + solvers. **Scaling/speed (landed):** the cost is the eigensolve, and oversampling + scales fine via ARPACK — the bottleneck was the dense-eigensolve threshold. + ``DENSE_EIG_MAX`` was lowered 256 → 50 (the measured dense/ARPACK crossover; dense + is O(N³), ARPACK ~flat in N on the banded laplacian), and the inner trust-region + solve was stopped from over-converging against the frozen field (tol 1e-6, + fewer field refreshes). line_PRA newton 47s → 14s, byte-identical. Optional + future refinements: an analytic coherent within-edge hole-burning integral (to + drop oversampling entirely) and a Jacobian-free amplitude update. ## Git / branch policy for this repo diff --git a/benchmark/bench_salt.py b/benchmark/bench_salt.py new file mode 100644 index 00000000..e046cc47 --- /dev/null +++ b/benchmark/bench_salt.py @@ -0,0 +1,235 @@ +"""Modal-intensity solver benchmark: linear vs full_salt_newton. + +Compares the two L--I (intensity-vs-pump) solvers (issue #42) on **speed** and +**accuracy**: + +* ``linear`` — the near-threshold SALT model (one pump-independent + competition matrix, a linear solve, piecewise-linear curves). +* ``full_salt_newton`` — the operator-level nonlinear SALT (solves the + saturated eigenproblem per pump; reduces to ``linear`` at threshold). + +The expensive, shared pipeline steps (passive modes, pump, trajectories, +thresholds, competition matrix) are run once through the normal cached pipeline; +only the modal-intensity step is swapped between solvers, so the table isolates +the cost and the drift of the two intensity models. + +Usage:: + + python benchmark/bench_salt.py # default: line_PRA + python benchmark/bench_salt.py examples/line_PRA/config.yaml + +Writes ``benchmark/bench_salt_ll.pdf`` (overlaid L--I curves) and +``benchmark/bench_salt_newton.pdf`` (the newton study). Requires the example to +be runnable from its own directory. +""" + +from __future__ import annotations + +import os +import sys +from pathlib import Path +from time import perf_counter + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from netsalt.config_loader import load_config +from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, +) +from netsalt.pipeline import ( + _attach_pump_to_graph, + step_compute_mode_competition_matrix, + step_compute_mode_trajectories, + step_create_pump_profile, + step_create_quantum_graph, + step_find_passive_modes, + step_find_threshold_modes, + step_scan_frequencies, +) + +HERE = Path(__file__).resolve().parent +REPO = HERE.parent + + +def time_block(): + class _T: + def __enter__(self): + self.t0 = perf_counter() + return self + + def __exit__(self, *a): + self.seconds = perf_counter() - self.t0 + + return _T() + + +def _ll_curve(df): + """Extract ``(pumps, per_mode, total)`` from an intensities dataframe.""" + cols = [c for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + pumps = np.array([c[1] for c in cols], dtype=float) + order = np.argsort(pumps) + pumps = pumps[order] + data = df[[cols[i] for i in order]].to_numpy(dtype=float) # (n_modes, n_pumps) + total = np.nansum(data, axis=0) + return pumps, data, total + + +def _shared_steps(p, lasing_modes_id): + """Run the cached pipeline up to (and including) the competition matrix.""" + qg = step_create_quantum_graph(p) + qualities = step_scan_frequencies(p, qg) + passive_modes_df = step_find_passive_modes(p, qg, qualities) + pump = step_create_pump_profile(p, qg, passive_modes_df, lasing_modes_id) + trajectories_df = step_compute_mode_trajectories(p, qg, passive_modes_df, pump, lasing_modes_id) + threshold_modes_df = step_find_threshold_modes(p, qg, trajectories_df, pump, lasing_modes_id) + competition = step_compute_mode_competition_matrix( + p, qg, threshold_modes_df, pump, lasing_modes_id + ) + qg = _attach_pump_to_graph(p, qg, pump) + return qg, threshold_modes_df, competition + + +def _solvers(p): + """Map method name -> zero-arg callable returning an intensities dataframe.""" + d0_max = p.get("intensities_D0_max") or p.get("D0_max", 0.1) + steps = p.get("salt_D0_steps", 30) + + def run_linear(qg, tdf, comp): + return compute_modal_intensities(tdf.copy(), d0_max, comp) + + def run_newton(qg, tdf, comp): + return compute_modal_intensities_full_salt_newton(qg, tdf.copy(), d0_max, D0_steps=steps) + + return {"linear": run_linear, "full_salt_newton": run_newton} + + +def benchmark(config_path: Path): + config_path = config_path.resolve() + workdir = config_path.parent + cwd = os.getcwd() + os.chdir(workdir) + try: + p = load_config(config_path.name) + lasing_modes_id = p.get("lasing_modes_id") + print(f"=== {config_path.relative_to(REPO)} ===") + + qg, tdf, comp = _shared_steps(p, lasing_modes_id) + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + n_lasing = int(np.sum(thresholds < np.inf)) + print(f"graph: {len(qg)} nodes, {len(qg.edges)} edges; {n_lasing} lasing modes") + + solvers = _solvers(p) + results = {} + print( + f"\n{'method':>16} | {'time (s)':>8} | {'n_active':>8} | {'tot@max':>10} | {'Δ vs lin':>9}" + ) + print("-" * 64) + + linear_total_max = None + for name in ("linear", "full_salt_newton"): + with time_block() as t: + df = solvers[name](qg, tdf, comp) + pumps, per_mode, total = _ll_curve(df) + results[name] = (pumps, per_mode, total) + tot_max = total[-1] if len(total) else float("nan") + n_active = int(np.sum(np.nan_to_num(per_mode[:, -1]) > 0)) if per_mode.size else 0 + if name == "linear": + linear_total_max = tot_max + delta = "-" + else: + delta = f"{(tot_max - linear_total_max):+.3e}" + print(f"{name:>16} | {t.seconds:>8.3f} | {n_active:>8} | {tot_max:>10.3e} | {delta:>9}") + + _plot_ll(results, HERE / "bench_salt_ll.pdf") + print(f"\nwrote {(HERE / 'bench_salt_ll.pdf').relative_to(REPO)}") + + _newton_study(qg, tdf, p, HERE / "bench_salt_newton.pdf") + finally: + os.chdir(cwd) + + +def _newton_study(qg, tdf, p, out, steps=8): + """Operator-level full-SALT Newton vs the linear model. + + Highlights two things: (1) the dominant mode's onset slope reduces to the + linear ``1/(T_μμ·D0_thr)`` near threshold, and (2) above threshold full SALT + deviates -- bent curves and competition-shifted secondary modes. With the + within-edge hole burning resolved (default auto-oversampling) the lasing set + agrees with the linear/competition-matrix count (e.g. both modes on + ``line_PRA``, matching Ge-Chong-Stone Eq. 28); a too-coarse mesh over-clamps + and spuriously drops modes. The Newton solve is expensive (a nested + frequency/profile + amplitude solve per pump, on an oversampled graph), so it + runs on a coarse ``steps`` grid. + """ + d0_max = p.get("intensities_D0_max") or p.get("D0_max", 0.1) + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + if not np.any(thresholds < np.inf): + return + target = int(np.argmin(thresholds)) + d0_thr = float(thresholds[target]) + + t_lin = compute_mode_competition_matrix(qg, tdf) + lin_df = compute_modal_intensities(tdf.copy(), d0_max, t_lin) + lin_last = np.nan_to_num(lin_df["modal_intensities"].to_numpy()[:, -1]) + lin_slope = 1.0 / (t_lin[target, target] * d0_thr) + + with time_block() as t: + df = compute_modal_intensities_full_salt_newton(qg, tdf.copy(), d0_max, D0_steps=steps) + sub = df["modal_intensities"] + pumps = np.array(sorted(sub.columns)) + new_last = np.nan_to_num(sub.to_numpy()[:, -1]) + a = sub.loc[target, pumps].to_numpy(dtype=float) + above = pumps > d0_thr + 1e-9 + newton_slope = a[above][0] / (pumps[above][0] - d0_thr) if above.any() else float("nan") + + print(f"\noperator-level Newton full SALT ({t.seconds:.0f} s, {len(pumps)} pumps):") + print(f" dominant-mode onset slope newton/linear: {newton_slope / lin_slope:.3f} (1.0 = ok)") + print(f" linear lases modes : {sorted(int(i) for i in np.where(lin_last > 0)[0])}") + print(f" newton lases modes : {sorted(int(i) for i in np.where(new_last > 1e-6)[0])}") + print( + " (counts should agree with resolved hole burning; full SALT bends the curves above threshold)" + ) + + plt.figure(figsize=(6, 4)) + plt.plot( + pumps, np.clip(lin_slope * (pumps / d0_thr - 1.0), 0, None), "--", label="linear (mode)" + ) + plt.plot(pumps, a, "o-", ms=3, label="full_salt_newton (mode)") + plt.xlabel("pump $D_0$") + plt.ylabel(f"dominant-mode intensity (mode {target})") + plt.legend() + plt.title("Operator-level Newton vs linear (dominant mode)") + plt.tight_layout() + plt.savefig(out) + plt.close() + print(f"wrote {out.relative_to(REPO)}") + + +def _plot_ll(results, out): + plt.figure(figsize=(6, 4)) + for name, (pumps, _per_mode, total) in results.items(): + plt.plot(pumps, total, marker="o", ms=3, label=name) + plt.xlabel("pump $D_0$") + plt.ylabel("total modal intensity") + plt.legend() + plt.title("L--I curves: linear vs full_salt_newton") + plt.tight_layout() + plt.savefig(out) + plt.close() + + +def main(argv=None): + argv = sys.argv[1:] if argv is None else argv + config = Path(argv[0]) if argv else REPO / "examples" / "line_PRA" / "config.yaml" + benchmark(config) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/doc/cf_states_design.md b/doc/cf_states_design.md new file mode 100644 index 00000000..8e1e302f --- /dev/null +++ b/doc/cf_states_design.md @@ -0,0 +1,130 @@ +# Constant-flux (CF) state SALT solver — design and validated foundation + +> **Correction (this PR): `full_salt_newton` already runs on the real buffon.** +> The premise below — that buffon needs a CF solver because oversampling makes a +> ~77 000-node operator — was **wrong**. `_auto_oversample_size` bounds the +> oversampled graph by `node_cap` (default 3000), so newton runs on the real +> buffon in **~12 s** (default cap) and **~6 min uncapped at λ/4** (≈30k nodes), +> lasing its co-lasing modes. The node count is `node_cap`, not +> `inner_length / (λ/12)`. The `oversample_node_cap` / `oversample_resolution` +> knobs now expose the speed/accuracy trade. So the CF route is **not needed for +> feasibility**; it is kept below only as a (validated) curiosity and a possible +> future accuracy/scaling refinement. The rest of this note is retained for the +> record, but its motivation no longer holds. + +Status: **CF foundation validated; not required — buffon works via bounded +oversampling.** + +## Why CF states (the original, now-corrected motivation) + +`full_salt_newton` resolves the within-edge hole burning by **oversampling** the +graph to ~λ/12 sub-edges — *bounded by `node_cap`*. The reasoning below assumed +the cap did not exist (an error); it is kept for the record. + +Two cheaper routes were tried and **rejected**, each for a concrete measured +reason: + +* **Analytic per-edge gain (oversample-free).** Replace oversampling with the + closed-form within-edge integral `⟨1/(1+Γa|E(x)|²)⟩` by Gauss-Legendre + quadrature of the exact plane-wave field. *Validated correct for a single + mode* (matches oversample on line_PRA in the single-mode regime, 5.8× faster). + But it **cannot discriminate co-lasing modes**: on line_PRA at D0 = 1.27 it + dumped all intensity on the dominant mode (0.30) and suppressed the second + (0.0) where oversample lases two (0.21 / 0.099). Folding the gain into one + *constant per edge* — even the exact edge-average — washes out the within-edge + standing-wave structure that lets two modes share the gain (each saturates it + where *its* field peaks). That structure lives in the operator; a per-edge + constant cannot carry it. *Series expansion of `1/(1+s)` also fails* — it + diverges (s > 1 on antinode edges even at moderate amplitude). +* **Coarser oversampling.** Discriminating modes that vary on scale λ needs the + gain resolved on scale λ *in the operator*, so the node count stays O(L/λ) + regardless of how exactly `|E|²` is computed. No real win. + +The CF basis is the principled fix: it represents the field in **exact graph +eigenfunctions** (plane waves per edge), so (i) overlap integrals are analytic +(the competition-matrix machinery), **no oversampling**, and (ii) the basis +carries the full within-edge spatial structure, so it **resolves multimode**. + +## The quantum-graph CF eigenproblem + +Ge-Chong-Stone (PRA 82, 063824) Eq. (17) defines the threshold CF (TCF) states +at a fixed real frequency k: + + [∇² + ε(r) k² + ηₙ(k) F(r) k²] uₙ(r, k) = 0, outgoing BC, + +where `F` is the pump profile and `ηₙ` (complex, `Im ηₙ < 0`) is the **CF +eigenvalue** — the scale of amplifying dielectric needed for that state to reach +threshold at k. On a quantum graph this means the per-edge wavenumber is + + k_edge = k · √(ε_e + ηₙ · pump_e), + +i.e. the secular matrix `L(k, η)` built with the per-edge dielectric +`ε + η·pump` (plain dispersion) is **singular** for the CF eigenvalues `ηₙ(k)`. +This is a *nonlinear eigenproblem in η* at fixed real k — the same shape as the +passive mode search (nonlinear in k), so the existing Beyn machinery applies in +the η-plane. + +### Validated foundation (line_PRA) + +The **threshold** CF eigenvalue must equal the threshold gain `γ(k_thr)·D0_thr`, +because at threshold the amplifying dielectric *is* the actual gain and the CF +state *is* the threshold lasing mode (the paper's Fig. 3). Measured on +line_PRA's dominant mode (`/tmp` prototype, reproducible): + + k_thr = 15.44376, D0_thr = 0.61069, γ(k_thr) = 0.14475 − 0.97859j + expected η = γ·D0_thr = 0.08840 − 0.59762j + scanned argmin |λ₁(η)| at η = 0.08840 − 0.59762j (|λ₁| = 6e-7) + agreement: 0.00% + +So the operator `construct_laplacian(k_thr, graph_with dielectric ε+η·pump)` is +singular *exactly* at `η = γ·D0_thr` — the CF formulation is correct and +computable on the original 11-node graph, no oversampling. + +## Plan (the multi-PR work ahead) + +1. **CF basis finder — the hard step.** Two approaches were tried and **both + fail**, which is the central open problem: + + * *Beyn in the η-plane* — run the `contour.py` Cauchy-moment extraction on + `L(η) = construct_laplacian(k, dielectric = ε + η·pump)` instead of `L(k)`. + Tested on line_PRA: it returns a wrong root (Im η > 0, absorbing) and + **misses the validated threshold root even with a tight contour around + it** (found 0). Diagnosis: on a graph the secular matrix depends on η + through `√(ε + η·pump)` in the per-edge wavenumber, so `L(η)` is **not + analytic** in η (branch points at `η = −ε_e/pump_e`) and the + near-null behaviour around a CF eigenvalue is *soft*, not a clean simple + pole — Beyn's contour integral, which assumes a meromorphic `L(η)⁻¹`, + does not extract it. (This is the key difference from the *continuous* CF + operator of Eq. 17, which is **linear** in η.) + * *Coarse grid scan of `|λ₁(η)|`* — the η-roots are isolated points; a 61×61 + box found only the threshold one (which was centred), none of the basis. + + So the basis finder is genuine research, not plumbing. Candidate routes: + (a) **nonlinear root-finding in η** (Newton / secant on `λ₁(η)` or + `det L(η)`) seeded from physical guesses (passive modes near k, UCF + estimates `ηₙ = c(Kₙ²/k² − 1)` from Eq. 21) — handles non-analyticity since + it never integrates around the root; (b) a **linear-in-η reformulation** — + the edge ODE `u'' + k²(ε+η pump)u = 0` *is* linear in η, so a spatial + discretisation gives a standard generalised eigenproblem + `(K − k²M_ε)u = k²η M_pump u` solving all CF states at once, but that + reintroduces a mesh (the cost we are trying to avoid); (c) a **hybrid** — + coarse mesh only to seed guesses for the exact nonlinear η-solve. + Route (a) is the most promising and is the next concrete experiment. +2. **Self-orthogonality + overlaps.** CF states obey the paper's self-orthogonality + relation (Eq. 20); the hole-burning overlaps `∫ uₙ uₘ* / (1 + Σ |Ψ|²)` over + each edge are products of plane waves → closed-form exponential integrals + (extend `_compute_mode_competition_element`). No oversampling. +3. **SALT in the CF basis (Eq. 28).** Expand each lasing mode `Ψμ = Σₙ aₙᵘ uₙ(kμ)`; + the SALT condition becomes the nonlinear fixed-point map `Tₙₙ'(k) a = a` + (Eq. 28) in the CF coefficients, with the gain-clamping nonlinearity in the + overlaps. Solve `(kμ, aₙᵘ)` self-consistently — the lasing map's fixed points, + no spatial discretization. +4. **Validation milestones.** (a) threshold CF = threshold mode (✓ done); + (b) CF basis reproduces a passive mode by expansion; (c) single-mode L–I + matches oversample on line_PRA; (d) **two-mode** L–I reproduces Ge Fig. 6 + (the test the per-edge-average failed); (e) run the real buffon at its native + node count. + +The payoff: steps 1–3 keep the operator at the original node count, so the real +buffon runs at 96 nodes instead of 77 000 — and unlike the per-edge-average +shortcut, the CF basis resolves co-lasing modes. diff --git a/doc/source/lasing.rst b/doc/source/lasing.rst index 38d7ccb6..d9799ec9 100644 --- a/doc/source/lasing.rst +++ b/doc/source/lasing.rst @@ -186,8 +186,8 @@ Approximations and validity across the package (loss = positive imaginary part); the threshold and competition expressions are consistent with it. -Relation to full SALT (future work) ------------------------------------ +Relation to full SALT +--------------------- The model above is the *linearized* (near-threshold) limit of SALT. The full SALT equation for each lasing mode :math:`\Psi_\mu` at real frequency @@ -241,8 +241,10 @@ improved independently. Relaxing them defines a hierarchy: :math:`\sum_\nu T_{\mu\nu} I_\nu = D_0/D_0^{\mathrm{thr}}_\mu - 1` - competition to first order; exact at threshold * - **Self-consistent linearized** (relax the frozen profile only) - - fixed-point loop: solve :math:`I` → rebuild :math:`T` from the profiles - at the *current* :math:`D_0` → repeat + - same event-driven sweep, but :math:`T` is rebuilt from the profiles at + each operating :math:`D_0` instead of held fixed at threshold. With linear + saturation :math:`T` depends only on the pump, so no inner fixed point is + needed. - profile deformation, gain guiding, frequency pulling; saturation still linear * - **Full SALT** (relax both) @@ -256,8 +258,164 @@ Quantum graphs are a favourable setting for the full solve: the per-edge field is two analytic plane waves, the secular matrix :math:`L(k) = B^{\mathsf T} W^{-1} B` is already assembled, and the saturated-gain integrals are the same closed-form edge integrals used in the competition -matrix (``_compute_mode_competition_element``). A full-SALT solver would fold the -hole-burning denominator into the gain term of -:func:`~netsalt.quantum_graph.construct_laplacian` and Newton-solve over the -modal amplitudes and frequencies at each pump, reusing those overlaps, with the -linearized model as the threshold-limit check. This is tracked as future work. +matrix (``_compute_mode_competition_element``). + +How many modes lase? Gain clamping vs. competition +-------------------------------------------------- + +The most visible difference between the solvers is **how many modes they lase**, +and it comes straight from the hole-burning denominator. Once mode :math:`\mu` +lases it **clamps** the saturated gain at its own threshold level. A second mode +:math:`\nu` keeps lasing only if it still has net gain *after* that clamping, +which depends on how much its intensity :math:`|E_\nu|^2` **overlaps** mode +:math:`\mu`'s spatial hole — exactly the off-diagonal :math:`T_{\mu\nu}` relative +to the self-saturation :math:`T_{\mu\mu}`: + +* **Strong overlap** (modes share the same region of the graph, e.g. a short + cavity with a narrow gain line) — the second mode is starved → *winner-take-all* + single-mode lasing. +* **Weak overlap** (modes occupy different regions / have distinct standing-wave + patterns, e.g. a broad gain line exciting well-separated modes) — the second + mode finds gain the first did not burn → *multimode* lasing. + +The solvers treat this clamping at different levels of fidelity: + +* ``linear`` switches a mode on when its fixed-:math:`T` *interacting threshold* is + crossed and never re-tests it. It captures the near-threshold competition exactly + (that is what :math:`T` is) but freezes the mode profiles, so above threshold it + misses how the deepening holes reshape the competition. +* ``full_salt_newton`` imposes the exact operator condition (the saturated operator + singular at real :math:`k` with :math:`a\ge 0`) — the *operator-level* SALT, rather + than the SPA matrix equation ``D0/D0_thr - 1 = Σ_ν Γ_ν χ_μν I_ν`` of + Ge–Chong–Stone. It contributes a self-consistent gain-clamping + **active set** and the lasing **frequencies** :math:`k_\mu` that the + competition-matrix solver cannot: on ``line_PRA`` it lases the **two** modes of + Ge–Chong–Stone (PRA 82, 063824, Eq. 28). Its amplitude is put in the linear unit + by an *analytic* onset scale + (first-order perturbation of the saturated operator's lasing condition, built + from the same coherent overlap integrals as ``pump_linear``), so it + **reduces to** ``linear`` at threshold -- the scale comes out ≈ 1, making the + near-threshold agreement a genuine prediction rather than a calibration. + +.. note:: + + **Validated against Ge–Chong–Stone Fig. 6.** Above threshold ``full_salt_newton`` + gives the genuine full-SALT correction beyond the SPA: when a second mode turns + on, the dominant mode picks up a **negative kink** and is suppressed *below* the + SPA (its gain is stolen), while the second mode sits *above* the SPA, the two + nearly cancelling in the total. On ``line_PRA`` the per-mode intensities track the + digitized exact-SALT curves of Fig. 6 to a few percent (dominant 0.21 vs 0.205, + second 0.10 vs 0.108 at :math:`D_0 = 1.27`), and the single-mode regime reduces to + the SPA. This requires *measuring* the onset slope: an earlier analytic-only scale + over-shot it by ~10–30 %, lifting the whole curve above the exact result. The + competition-matrix solvers remain the cheaper first pass; ``full_salt_newton`` + adds the operator-level above-threshold correction. + +.. warning:: + + The operator-level hole burning samples :math:`|E_\nu(x)|^2` per edge. With the + bare edges (one sample per edge) the per-edge **mean** over-estimates the mode + overlap — it washes out the standing-wave nodes/antinodes where the coherent + competition is weak — and **over-clamps**, spuriously suppressing co-lasing + modes. On ``line_PRA`` this made ``full_salt_newton`` lase one mode where Ge + Eq. 28 and the competition matrix lase two. ``oversample_size=None`` therefore + auto-picks a wavelength-resolving sub-edge size + (:func:`~netsalt.modes._auto_oversample_size`); resolving the standing wave + removes the over-clamping and recovers the correct count. The oversampled graph + is larger, but the eigensolve stays cheap because it runs through ARPACK + shift-invert (which is ~flat in the node count on the banded quantum-graph + laplacian; see ``DENSE_EIG_MAX``), and the oversampling is bounded by + ``oversample_node_cap`` (default 3000), so the operator stays a few thousand + nodes regardless of the cavity length -- the **real buffon network runs in + ~12 s** (default cap) and ~6 min uncapped at ``λ/4``. Raise + ``oversample_node_cap`` / ``oversample_resolution`` for within-edge accuracy + on large graphs at higher cost. ``linear`` remains the cheaper first pass for + the lasing count. + +.. note:: + + Multimode lasing is demonstrated in + ``examples/two_ring``: two **detuned** rings + (different sizes) joined by a bridge. The detuning localises each mode onto + one ring -- identical rings would give symmetric/antisymmetric modes spread + over both, with high overlap -- so with a narrow gain the modes barely + compete and ``full_salt_newton`` lases several at once (spread across the two + rings). ``examples/chaotic_ring`` shows + the same effect on a *single* small graph: one 14-node ring with six random + chords (the buffon mechanism shrunk down). The chords close extra loops, so + the spectrum is dense and the modes localise on different loops; with a + narrow gain on a four-mode cluster ``full_salt_newton`` lases four. Getting + there did need a tight ``k``-window (a loose one let the trust region collapse + the multimode set to one mode by drifting a mode's ``k`` to a spurious + ``a = 0`` root). Multimode remains the more delicate path, so treat it as + experimental and sanity-check the mode count. + +Selecting a solver +^^^^^^^^^^^^^^^^^^ + +Two solvers are available and chosen with the ``intensity_method`` config +key (default ``"linear"``), dispatched by +:func:`~netsalt.pipeline.step_compute_modal_intensities`. (Two intermediate +solvers, ``self_consistent`` and ``full_salt``, which rebuilt or saturated the +competition matrix at the operating pump, were removed: their per-pump matrix +rebuild was ill-conditioned in exactly the strongly-multimode regime where they +would have added value over ``linear`` — non-deterministic run-to-run, modes +locking to equal intensities or collapsing to zero — and on weakly-competing +graphs they only track ``linear``.) + +``"linear"`` + :func:`~netsalt.modes.compute_modal_intensities` — the + near-threshold model described above. This is the default; the newton + solver reduces to it at threshold. Fast. +``"full_salt_newton"`` + :func:`~netsalt.modes.compute_modal_intensities_full_salt_newton` — + *experimental, operator-level.* Rather than saturating the competition + matrix, it solves the real nonlinear SALT eigenproblem: at each pump it finds, + for every lasing mode, ``(k_μ, a_μ)`` so the shared saturated operator + ``L_sat`` (:func:`~netsalt.physics.dispersion_relation_pump_saturated`) is + singular at each real ``k_μ``. Two ingredients make it robust: + + * **Frozen-field trust-region solve** (:func:`~netsalt.modes._solve_active_set`) + -- for a fixed active set the saturated background fields are frozen while a + bounded trust-region least-squares solves all ``(k_μ, a_μ)``; the fields are + then refreshed and the step repeated. Freezing the field makes each residual a + single clean eigensolve (no inner fixed point), so the Jacobian is noise-free + -- the earlier decoupled solve, whose residual re-ran an inner fixed point, + chattered for several co-lasing modes. + * **Self-consistent active-set continuation** -- the pump is stepped up; the + confirmed lasing set is solved, modes whose amplitude vanishes are dropped, + and a candidate is added when it has net gain (``α < 0``) on the current + saturated background. The active set is found from the saturated operator, not + borrowed from the linear model. This is only faithful when the within-edge + hole burning is **resolved**: ``oversample_size=None`` auto-picks a + wavelength-resolving sub-edge size, without which the per-edge mean + over-clamps and spuriously drops co-lasing modes (see the warning above). + + Its amplitude is put in the **linear modal-intensity unit** by an *analytic* + onset scale (first-order perturbation of the saturated operator's lasing + condition; it comes out ≈ 1), so it **reduces to linear at + threshold** and agrees on the lasing count -- validated against Ge–Chong–Stone + (PRA 82, 063824) on ``line_PRA`` (both lase two modes). Above threshold it is + the *exact-spatial* SALT: the + dominant mode gets a negative kink (suppressed below the SPA when the second mode + steals gain) and the second mode sits above the SPA, tracking the **exact-SALT + data of Fig. 6 to a few percent**. The linear + solver remains the cheaper first pass. It is deterministic, path-independent, + never raises, and *expensive* (a nested per-pump solve on the oversampled graph), + so use a modest + ``salt_D0_steps``. + +``benchmark/bench_salt.py`` compares the two solvers on speed and accuracy: it +runs the shared pipeline once, swaps only the intensity step, writes overlaid +L–I curves, and contrasts the operator-level Newton solver with the linear model +(onset slope + which modes lase). The newton solver is exact at threshold, so +the linear model remains the threshold-limit check. + +The script-based example folders (one per graph) are self-contained worked +examples: the simple open cavities (``examples/line_fabry_perot``, +``examples/ring_leads``, ``examples/tree``) overlay the two methods' L–I curves +with a per-mode breakdown that makes the above-threshold bend-over explicit, and +the multimode graphs (``examples/two_ring``, ``examples/chaotic_ring``, +``examples/dense_ring``) probe the operator-level mode competition. Each +folder's ``run.py`` regenerates its figures (figures are not committed). + diff --git a/doc/source/physics.rst b/doc/source/physics.rst index 571a7e0a..bad98ee2 100644 --- a/doc/source/physics.rst +++ b/doc/source/physics.rst @@ -12,7 +12,14 @@ through the dielectric law (:func:`~netsalt.physics.dispersion_relation_dielectric`, :math:`k = \omega\sqrt{\varepsilon}/c`) to the pumped law (:func:`~netsalt.physics.dispersion_relation_pump`), which adds the gain term -:math:`\gamma(k)\,D_0\,\delta_\mathrm{pump}` under the square root. The gain +:math:`\gamma(k)\,D_0\,\delta_\mathrm{pump}` under the square root. The +nonlinear-SALT solvers use +:func:`~netsalt.physics.dispersion_relation_pump_saturated`, the same law with +that gain term replaced by a per-edge *saturated* effective pump +:math:`D_0^\mathrm{eff} = D_0\,\delta_\mathrm{pump} / (1 + \sum_\nu \Gamma_\nu +a_\nu |\Psi_\nu|^2)` carrying the spatial-hole-burning denominator (see +:doc:`lasing`); with the denominator set to one it reduces to the pumped law. +The gain itself is the Lorentzian :func:`~netsalt.physics.gamma`, centred on the atomic transition wavenumber ``k_a`` with linewidth ``gamma_perp``; ``D0`` is the pump strength and ``params["pump"]`` the per-edge pump profile. diff --git a/examples/README b/examples/README index 260772a5..98354379 100644 --- a/examples/README +++ b/examples/README @@ -1,9 +1,18 @@ -This folder contains examples of netSALT runs found in the accompanying paper. +This folder contains examples of netSALT runs found in the accompanying paper, +plus script-based per-graph examples of the modal-intensity solvers. In each, run ``` ./run.sh ``` -to run the workflow, and adjust parameters in `config.yaml` for the simulation -parameters. Most variants inherit a shared `_base.yaml` via `defaults:` — -override only the keys that differ. +For the config-driven examples (buffon, ring, wheel, directed, line_PRA, +transfer) this runs the pipeline workflow; adjust parameters in `config.yaml` +(most variants inherit a shared `_base.yaml` via `defaults:` — override only +the keys that differ). + +The script-based examples (`line_fabry_perot`, `ring_leads`, `tree`, +`two_ring`, `chaotic_ring`, `dense_ring`; shared helpers in `_common.py`) +build one graph each in memory, compare the `linear` and `full_salt_newton` +intensity solvers, and write their figures into their own folder. Figures are +not committed — re-run the script to reproduce them. See +`doc/source/lasing.rst` for the physics. diff --git a/examples/_common.py b/examples/_common.py new file mode 100644 index 00000000..9d0194a7 --- /dev/null +++ b/examples/_common.py @@ -0,0 +1,338 @@ +"""Shared machinery for the script-based per-graph examples. + +The per-graph example folders (``line_fabry_perot/``, ``ring_leads/``, +``tree/``, ``two_ring/``, ``chaotic_ring/``, ``dense_ring/``) each hold a +``run.py`` that builds one graph in memory, runs the two modal-intensity +solvers (``linear`` and ``full_salt_newton``) and writes its figures *into +that folder* (figures are gitignored -- re-run the script to reproduce them). +This module carries what they share: the open-cavity parameter set, the +pipeline runner, curve extraction, and the standard geometry/per-mode/total +comparison figure used by the three simple cavities. +""" + +from __future__ import annotations + +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np + +import netsalt +from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + find_passive_modes, + find_threshold_lasing_modes, + pump_trajectories, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length + +# Shared physics / search settings for the three simple open cavities (a gain +# line centred at ``k_a`` and a scan window straddling it). Kept small so each +# run.py finishes in about a minute on one core. +OPEN_CAVITY_PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 15.0, + "gamma_perp": 3.0, + "k_min": 12.0, + "k_max": 18.0, + "k_n": 100, + "alpha_min": 0.0, + "alpha_max": 1.0, + "alpha_n": 30, + "quality_threshold": 1.0e-4, + "search_stepsize": 0.01, + "max_steps": 10000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "n_workers": 1, + "D0_max": 1.4, + "D0_steps": 10, + "salt_D0_steps": 10, + "dielectric_params": { + "method": "uniform", + "inner_value": 9.0, # epsilon = n^2, n = 3 inside the cavity + "outer_value": 1.0, # the leads are vacuum + "loss": 0.0, + }, +} +D0_MAX = 1.4 # max pump for the simple-cavity L--I sweeps + + +def quantum_graph(nx_graph, positions, total_length, **overrides): + """Wrap a networkx graph as a pumped, open netSALT quantum graph. + + ``overrides`` patch the shared :data:`OPEN_CAVITY_PARAMS`. + """ + g = nx.convert_node_labels_to_integers(nx_graph) + params = dict(OPEN_CAVITY_PARAMS) + params.update(overrides) + create_quantum_graph(g, params, positions=positions) + set_total_length(g, total_length) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def threshold_modes(graph, method="grid", pump=None): + """Shared pipeline: passive modes -> pump -> trajectories -> thresholds. + + ``method="grid"`` scans the frequency grid first (the simple cavities); + ``method="contour"`` uses Beyn's contour search (the ring graphs). + ``pump`` is the per-edge pump profile (default: uniform on the inner + edges); trajectories and thresholds are computed under it. + """ + if method == "grid": + qualities = netsalt.scan_frequencies(graph) + passive = netsalt.find_passive_modes( + graph, qualities, method="grid", min_distance=2, threshold_abs=0.1 + ) + else: + passive = find_passive_modes(graph, method="contour") + if len(passive) == 0: + raise SystemExit("no passive modes found") + if pump is None: + pump = np.array([1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()]) + graph.graph["params"]["pump"] = np.asarray(pump, dtype=float) + trajectories = pump_trajectories(passive, graph, return_approx=True) + return find_threshold_lasing_modes(trajectories, graph) + + +def curves(df): + """(pumps, intensities[n_modes, n_pumps]) from an L--I dataframe.""" + cols = np.array( + sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + ) + data = np.nan_to_num(df[[("modal_intensities", c) for c in cols]].to_numpy(dtype=float)) + return cols, data + + +def endpoint(df): + """Modal intensities at the largest pump of an L--I dataframe.""" + return curves(df)[1][:, -1] + + +def linear_on_grid(tdf, competition, grid, n_modes): + """Sample the (exact, piecewise-linear) linear model on a uniform pump grid. + + The event-driven sweep only stores points at activation events; endpoint- + sampling a uniform grid simply draws the same line at the plot resolution. + """ + out = np.zeros((n_modes, grid.size)) + for j, d0 in enumerate(grid): + out[:, j] = endpoint(compute_modal_intensities(tdf.copy(), d0, competition)) + return out + + +def draw_geometry(ax, graph, name): + """Cavity (inner) edges grey, leads red-dashed, using the stored positions.""" + pos = {n: np.asarray(graph.nodes[n]["position"], dtype=float) for n in graph.nodes} + for u, v in graph.edges(): + x = [pos[u][0], pos[v][0]] + y = [pos[u][1], pos[v][1]] + if graph[u][v]["inner"]: + ax.plot(x, y, color="0.55", lw=2.2, zorder=1) + else: + ax.plot(x, y, color="crimson", lw=2.0, ls="--", zorder=1) + for n in graph.nodes: + ax.scatter(*pos[n], s=80, color="white", edgecolor="black", zorder=3) + ax.set_aspect("equal") + ax.axis("off") + ax.set_title(f"{name} ({len(graph)} nodes)") + + +def compare_and_plot( + graph, name, outdir, d0_max=D0_MAX, d0_steps=26, passive_method="grid", pump=None, prefix="" +): + """Run linear + newton on ``graph`` and write the standard 3-panel figure. + + Panels: geometry | per-mode L--I (linear dashed, newton solid, colours keyed + by mode) | total L--I for both solvers. Writes ``/li_curves.png`` + and prints the lasing counts. Both solvers share the linear modal-intensity + unit (the newton amplitude reduces to it analytically at threshold), so the + curves overlay directly. + """ + import time + + outdir = Path(outdir) + tdf = threshold_modes(graph, method=passive_method, pump=pump) + n_modes = len(tdf) + competition = compute_mode_competition_matrix(graph, tdf) + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + first = float(thr[thr < np.inf].min()) + grid = np.linspace(first, d0_max, d0_steps) + t0 = time.perf_counter() + linear = linear_on_grid(tdf, competition, grid, n_modes) + t_linear = time.perf_counter() - t0 + t0 = time.perf_counter() + n_cols, newton = curves( + compute_modal_intensities_full_salt_newton(graph, tdf.copy(), d0_max, D0_steps=d0_steps) + ) + t_newton = time.perf_counter() - t0 + print(f"{name}: linear sweep {t_linear:.1f}s, newton sweep {t_newton:.1f}s") + + peak = max(linear.max(), newton.max(), 1e-9) + active = [m for m in range(n_modes) if max(linear[m].max(), newton[m].max()) > 1e-2 * peak] + n_lin = int(np.sum(linear[:, -1] > 1e-2 * peak)) + n_nwt = int(np.sum(newton[:, -1] > 1e-2 * peak)) + finite = np.sort(thr[thr < np.inf])[:6] + print(f"{name}: lowest thresholds {np.round(finite, 3)}") + print(f"{name}: linear lases {n_lin}, full_salt_newton lases {n_nwt}") + + cmap = plt.get_cmap("tab10") + fig, axes = plt.subplots(1, 3, figsize=(15, 4.2)) + draw_geometry(axes[0], graph, name) + for m in active: + col = cmap(m % 10) + axes[1].plot(grid, linear[m], "--", color=col, lw=1.3, alpha=0.8) + axes[1].plot(n_cols, newton[m], ".-", color=col, lw=1.8, ms=4, label=f"mode {m}") + axes[1].set_xlabel("pump $D_0$") + axes[1].set_ylabel("modal intensity") + axes[1].set_title(f"per mode: dashed = linear ({n_lin}), solid = newton ({n_nwt})") + if active: + axes[1].legend(fontsize=8, ncol=2) + axes[2].plot(grid, linear.sum(axis=0), "o-", ms=3, color="tab:blue", label="linear") + axes[2].plot(n_cols, newton.sum(axis=0), "o-", ms=3, color="tab:red", label="full_salt_newton") + axes[2].set_xlabel("pump $D_0$") + axes[2].set_ylabel("total modal intensity") + axes[2].set_title("total") + axes[2].legend(fontsize=8) + fig.suptitle(f"{name}: full_salt_newton vs linear", y=1.02) + fig.tight_layout() + out = outdir / f"{prefix}li_curves.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") + + ids_linear = [m for m in range(n_modes) if linear[m, -1] > 1e-2 * peak] + ids_newton = [m for m in range(n_modes) if newton[m, -1] > 1e-2 * peak] + mode_profile_figure( + graph, + tdf, + d0_max, + ids_linear, + ids_newton, + name, + outdir, + a0={m: newton[m, -1] for m in ids_newton}, + filename=f"{prefix}mode_profiles.png", + ) + return {"tdf": tdf, "linear": linear, "newton": newton, "grid": grid, "n_cols": n_cols} + + +def _draw_profile(ax, work, values, cmap, vmin, vmax): + """Edges of ``work`` coloured by the per-edge ``values`` (intensity map).""" + from matplotlib.collections import LineCollection + + pos = {n: np.asarray(work.nodes[n]["position"], dtype=float) for n in work.nodes} + segs = [(pos[u], pos[v]) for u, v in work.edges] + lc = LineCollection(segs, cmap=cmap, linewidths=3.0) + lc.set_array(np.asarray(values, dtype=float)) + lc.set_clim(vmin, vmax) + ax.add_collection(lc) + ax.autoscale() + ax.set_aspect("equal") + ax.axis("off") + return lc + + +def mode_profile_figure( + graph, tdf, d0_max, ids_linear, ids_newton, name, outdir, a0=None, filename="mode_profiles.png" +): + """Per-mode profile comparison: linear vs saturated (newton) vs difference. + + For every mode lasing under *either* solver this draws three columns: + + * **linear** -- the mode's intensity profile ``|E(x)|^2`` frozen at its own + threshold (the field the competition matrix is built from); + * **full_salt_newton** -- the profile of the same mode at the operating pump + ``d0_max`` under the shared *saturated* operator (all modes' spatial hole + burning included), obtained by re-solving the coupled ``(k, a)`` + equilibrium with the union of both lasing sets active, warm-started from + the linear profiles and the newton endpoint amplitudes; + * **difference** (saturated - linear) -- where the deepening holes reshape + the mode. This is the quantity the linear model freezes, so it is the + *reason* the two solvers can disagree on counts and interacting + thresholds: a suppressed/late mode is one whose saturated profile is + pushed off the gain, an early one finds gain the linear overlap + over-counted. + + All profiles share the ``int_pump |E|^2 = 1`` normalisation (per-edge values + on the oversampled work graph), so the columns are directly comparable. + Writes ``/mode_profiles.png``. + """ + import netsalt.modes as _m + import netsalt.quantum_graph as _qg + from netsalt.quantum_graph import graph_with_pump, oversample_graph + from netsalt.utils import from_complex + + ids = sorted({int(i) for i in ids_linear} | {int(i) for i in ids_newton}) + if not ids: + print(f"{name}: no lasing modes, skipping mode-profile figure") + return + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + tms = tdf["threshold_lasing_modes"].to_numpy() + + size = _m._auto_oversample_size(graph, tdf) + work = oversample_graph(graph, size) if size else graph + pump = np.asarray(work.graph["params"]["pump"], dtype=float) + pump_mask = _m._get_mask_matrices(work.graph["params"])[1] + + saved = _qg.DENSE_EIG_MAX + _qg.DENSE_EIG_MAX = min(saved, _m.NEWTON_DENSE_EIG_MAX) + try: + modes0, fields0 = [], [] + for i in ids: + mode = np.asarray(from_complex(tms[i]), dtype=float) + fields0.append( + _m._single_mode_field_intensity( + graph_with_pump(work, float(thr[i])), mode, pump_mask + ) + ) + modes0.append(mode) + start = ( + [max(float(a0.get(i, 1.0)), 1e-3) for i in ids] if a0 is not None else [1.0] * len(ids) + ) + modes, fields, amps, _ = _m._solve_active_set( + work, modes0, fields0, start, float(d0_max), pump, pump_mask, 30, 42 + ) + finally: + _qg.DENSE_EIG_MAX = saved + + fig, axes = plt.subplots(len(ids), 3, figsize=(13, 3.4 * len(ids)), squeeze=False) + for j, i in enumerate(ids): + row = j + lin, sat = np.asarray(fields0[j]), np.asarray(fields[j]) + vmax = max(lin.max(), sat.max()) + diff = sat - lin + dmax = max(abs(diff).max(), 1e-12) + lases = ("linear" if i in ids_linear else "") + ( + ("+newton" if i in ids_linear else "newton") if i in ids_newton else "" + ) + lc0 = _draw_profile(axes[row, 0], work, lin, "viridis", 0.0, vmax) + axes[row, 0].set_title( + f"mode {i} (k={modes0[j][0]:.3f}, thr={thr[i]:.3g}, lases: {lases})\n" + f"linear profile (at own threshold)", + fontsize=9, + ) + lc1 = _draw_profile(axes[row, 1], work, sat, "viridis", 0.0, vmax) + axes[row, 1].set_title( + f"saturated profile at D0={d0_max:g} (newton, a={amps[j]:.3g})", fontsize=9 + ) + lc2 = _draw_profile(axes[row, 2], work, diff, "coolwarm", -dmax, dmax) + axes[row, 2].set_title("difference (saturated - linear)", fontsize=9) + for lc, ax in ((lc0, axes[row, 0]), (lc1, axes[row, 1]), (lc2, axes[row, 2])): + fig.colorbar(lc, ax=ax, fraction=0.045, pad=0.02) + fig.suptitle(f"{name}: hole burning reshapes the mode profiles", y=1.0) + fig.tight_layout() + out = Path(outdir) / filename + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") diff --git a/examples/buffon/buffon_multimode/_base.yaml b/examples/buffon/buffon_multimode/_base.yaml index 93dc2842..38728949 100644 --- a/examples/buffon/buffon_multimode/_base.yaml +++ b/examples/buffon/buffon_multimode/_base.yaml @@ -4,6 +4,11 @@ n_workers: 10 D0_max: 0.02 intensities_D0_max: 0.02 +# Modal-intensity solver: "linear" (default, fast) or "full_salt_newton" +# (operator-level SALT; expensive on large multimode graphs). +# See doc/source/lasing.rst. +intensity_method: linear + pump_mode: optimized plot_threshold_mode_ids: [1, 5, 20] diff --git a/examples/buffon/buffon_multimode/buffon_2modes_newton/README.md b/examples/buffon/buffon_multimode/buffon_2modes_newton/README.md new file mode 100644 index 00000000..532dd3c6 --- /dev/null +++ b/examples/buffon/buffon_multimode/buffon_2modes_newton/README.md @@ -0,0 +1,16 @@ +# Buffon network with the operator-level solver + +Same pump-optimised two-mode buffon as `../buffon_2modes`, but with +`intensity_method: full_salt_newton` — operator-level full SALT on the **real +buffon network**, not a shrunk stand-in. + +This works because `full_salt_newton` bounds its oversampled operator by +`intensity_oversample_node_cap` (default 3000), so the eigensolve stays a few +thousand nodes *regardless of the cavity's physical length*. On the real buffon +it runs in ~12 s at the default cap and ~6 min at the raised cap set here +(`30000`, ≈ λ/4 resolution — the within-edge floor measured on `line_PRA` for +resolving co-lasing modes). Raise `intensity_oversample_node_cap` / +`intensity_oversample_resolution` for more accuracy at more cost; `linear` +remains the cheap first pass for the lasing count. + +Run with `bash run.sh` (needs `../../buffon_uniform/out` populated first). diff --git a/examples/buffon/buffon_multimode/buffon_2modes_newton/config.yaml b/examples/buffon/buffon_multimode/buffon_2modes_newton/config.yaml new file mode 100644 index 00000000..0f75883b --- /dev/null +++ b/examples/buffon/buffon_multimode/buffon_2modes_newton/config.yaml @@ -0,0 +1,15 @@ +defaults: ../_base.yaml + +# Same two-mode buffon as ../buffon_2modes (pump optimised to lase modes 1, 5), +# but solved with the operator-level full_salt_newton instead of linear -- the +# real buffon network runs through it in a few minutes (the oversampled +# operator is bounded by intensity_oversample_node_cap, NOT by the cavity +# length). Raise the cap / resolution for more within-edge accuracy at more cost. +plot_passive_edge_size: 100.0 +plot_passive_n_modes: 12 + +lasing_modes_id: [1, 5] + +intensity_method: full_salt_newton +salt_D0_steps: 10 +intensity_oversample_node_cap: 30000 # ~lambda/4 on this graph; default 3000 is faster/coarser diff --git a/examples/buffon/buffon_multimode/buffon_2modes_newton/run.sh b/examples/buffon/buffon_multimode/buffon_2modes_newton/run.sh new file mode 100755 index 00000000..617b308c --- /dev/null +++ b/examples/buffon/buffon_multimode/buffon_2modes_newton/run.sh @@ -0,0 +1,13 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +mkdir -p out figures + +# reuse the passive modes / graph computed by buffon_uniform +cp ../../buffon_uniform/out/passive_modes.h5 out/passive_modes.h5 +cp ../../buffon_uniform/out/qualities.h5 out/qualities.h5 +cp ../../buffon_uniform/out/quantum_graph.json out/quantum_graph.json + +python -m netsalt passive config.yaml +python -m netsalt lasing config.yaml diff --git a/examples/chaotic_ring/README.md b/examples/chaotic_ring/README.md new file mode 100644 index 00000000..7fe4e147 --- /dev/null +++ b/examples/chaotic_ring/README.md @@ -0,0 +1,13 @@ +# Single ring + random chords — multimode on one small graph + +A 14-node ring with 6 hard-coded random chords and two leads: the buffon +mechanism shrunk down. Each chord closes a loop, so the cavity supports many +interfering path lengths — a dense, irregular spectrum of spatially-distinct +modes (see the participation ratios the script prints). They burn their holes +in different places and co-lase: with a narrow gain on a four-mode cluster +`full_salt_newton` lases 4 modes where the clamping-free `linear` model lases +3. The figure shows the geometry, the full-range L–I, and a zoom on the onset +where the modes switch on in turn. + +`bash run.sh` writes `chaotic_ring_multimode.png` here. Figures are not +committed; re-run to reproduce. See `doc/source/lasing.rst` for the physics. diff --git a/examples/chaotic_ring/run.py b/examples/chaotic_ring/run.py new file mode 100644 index 00000000..b0764c6a --- /dev/null +++ b/examples/chaotic_ring/run.py @@ -0,0 +1,296 @@ +"""Genuine multimode lasing on a *single* small ring with random chords. + +``two_ring_multimode.py`` gets several modes to co-lase by joining two detuned +rings: the detuning localises each mode onto one ring, so they barely overlap. +But you do not need a multi-component graph for that -- a **single** ring with a +handful of chords (extra edges across it) already does the job, and on a much +smaller graph. This is the mechanism behind the multimode buffon networks, shrunk +to 14 nodes. + +Why chords give multimode lasing: + +* Each chord closes a new loop, so the cavity now supports many interfering + path-length combinations -> a **dense, irregular spectrum** (no clean + longitudinal comb). +* Those modes are **spatially distinct**: each one concentrates on a different + subset of loops/chords (see the per-mode *participation ratio* printed below), + so they burn their spatial holes in different places and leave gain for one + another -> they co-lase instead of competing winner-take-all. + +With a narrow gain centred on a cluster of these modes, ``full_salt_newton`` lases +**four** modes here. The chord layout is fixed (hard-coded, so the result is +reproducible regardless of the NumPy RNG); it came from a small seed scan picking +a graph whose spectrum has a clean four-mode cluster. + +**Which solvers to trust here.** This deep-multimode regime (4--5 strongly +clustered thresholds) is exactly where the cheap solvers part ways: + +* ``linear`` -- exact near-threshold model, gives clean piecewise-linear L--I + curves, but has **no gain clamping** (it never asks whether a lasing mode still + has net gain once the others saturate it), so its count can err either way -- + here it lases *one fewer* than newton, because its frozen-threshold competition + matrix over-estimates how strongly the cluster suppresses the fourth mode. +* ``full_salt_newton`` -- the operator-level solve stays smooth and physical and + imposes the exact self-consistent gain clamping. + +So this script plots ``linear`` (dashed) and ``full_salt_newton`` (solid) in +three panels: the graph geometry, the full-range L--I, and a **zoom on the +onset** (the four thresholds sit in ``0.016--0.026``, marked by dotted lines). +The zoom shows the modes switching on in turn, and in particular that newton +turns the fourth mode on *well above* its bare threshold -- gain clamping delays +it until enough pump is present, which is exactly what ``linear`` (no clamping) +misses. + +Run:: + + OMP_NUM_THREADS=1 python run.py + +Modes are found by Beyn's contour method (robust on this hand-built graph). +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np +from matplotlib.lines import Line2D + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import mode_profile_figure + +import netsalt +from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + find_passive_modes, + find_threshold_lasing_modes, + mode_on_nodes, + pump_trajectories, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import from_complex + +HERE = Path(__file__).resolve().parent +M = 14 # ring nodes +# Chords across the ring, plus a lead on node 0 and node M // 2 for output +# coupling. Hard-coded (from a seed scan) so the spectrum is reproducible. +CHORDS = [(0, 7), (1, 5), (2, 12), (3, 5), (3, 10), (7, 11)] + +# A narrow gain centred on the four-mode cluster near k = 2.81. Mode finding uses +# the contour method, so only k_min/k_max/alpha bounds matter for it. +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 2.81, + "gamma_perp": 0.10, + "k_min": 2.55, + "k_max": 3.25, + "alpha_min": -0.05, + "alpha_max": 0.25, + "n_workers": 1, + "n_modes_max": 40, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 0.5, + "D0_steps": 14, + "dielectric_params": {"method": "uniform", "inner_value": 9.0, "outer_value": 1.0, "loss": 0.0}, +} +D0_MAX = 0.5 +D0_STEPS = 22 # pump points over the full range +D0_ZOOM = 0.08 # onset window (the four thresholds sit in 0.016--0.026) +D0_STEPS_ZOOM = 22 # pump points within the zoom (fine, to resolve each turn-on) + + +def build_chaotic_ring(): + """A single ring with random chords and two leads (open cavity).""" + g = nx.cycle_graph(M) + g.add_edges_from(CHORDS) + g.add_edge(0, M) # lead + g.add_edge(M // 2, M + 1) # lead + pos = {i: [np.cos(2 * np.pi * i / M), np.sin(2 * np.pi * i / M)] for i in range(M)} + pos[M] = [1.6, 0.0] + pos[M + 1] = [-1.6, 0.0] + positions = np.array([pos[i] for i in range(len(g))]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 12.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def _participation(mode, graph): + """Effective number of nodes the mode lives on (1 / sum p_i^2).""" + weight = np.abs(mode_on_nodes(from_complex(mode), graph, check_quality=False)) ** 2 + prob = weight / weight.sum() + return 1.0 / np.sum(prob**2) + + +def _draw_geometry(ax, graph): + """Draw the ring/chord/lead geometry in the plane.""" + pos = {n: graph.nodes[n]["position"] for n in graph.nodes} + chord_set = {tuple(sorted(c)) for c in CHORDS} + for u, v in graph.edges(): + x = [pos[u][0], pos[v][0]] + y = [pos[u][1], pos[v][1]] + su, sv = min(u, v), max(u, v) + if sv >= M: # lead edge (nodes M, M + 1) + ax.plot(x, y, color="crimson", lw=2.0, ls="--", zorder=1) + elif (su, sv) in chord_set: + ax.plot(x, y, color="royalblue", lw=2.2, zorder=2) + else: # ring edge + ax.plot(x, y, color="0.6", lw=2.2, zorder=1) + for n in graph.nodes: + ax.scatter(*pos[n], s=220 if n < M else 160, color="white", edgecolor="black", zorder=3) + ax.text(pos[n][0], pos[n][1], str(n), ha="center", va="center", fontsize=7, zorder=4) + ax.legend( + handles=[ + Line2D([0], [0], color="0.6", lw=2.2, label="ring edge"), + Line2D([0], [0], color="royalblue", lw=2.2, label="random chord"), + Line2D([0], [0], color="crimson", lw=2.0, ls="--", label="output lead"), + ], + loc="upper right", + fontsize=7, + ) + ax.set_aspect("equal") + ax.axis("off") + ax.set_title(f"{M}-node ring + {len(CHORDS)} chords + 2 leads") + + +def _endpoint(df): + """Modal intensities at the largest pump of an L--I dataframe.""" + cols = sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + return np.nan_to_num(df[("modal_intensities", cols[-1])].to_numpy(dtype=float)) + + +def _curves(df): + """(pumps, intensities[n_modes, n_pumps]) from an L--I dataframe.""" + cols = np.array( + sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + ) + data = np.nan_to_num(df[[("modal_intensities", c) for c in cols]].to_numpy(dtype=float)) + return cols, data + + +def _linear_on_grid(tdf, competition, grid, n_modes): + """Sample the (exact, piecewise-linear) linear model on a uniform pump grid.""" + out = np.zeros((n_modes, grid.size)) + for j, d0 in enumerate(grid): + out[:, j] = _endpoint(compute_modal_intensities(tdf.copy(), d0, competition)) + return out + + +def _plot_li(ax, grid_lin, linear, n_cols, newton, thr, cmap, xmax=None): + """Overlay linear (dashed) and newton (solid) L--I curves, coloured by mode id.""" + n_modes = newton.shape[0] + active = [m for m in range(n_modes) if max(linear[m].max(), newton[m].max()) > 1e-3] + for m in active: + col = cmap(m % 10) + if thr[m] < (xmax if xmax is not None else np.inf): + ax.axvline(thr[m], color=col, ls=":", lw=0.8, alpha=0.6) + ax.plot(grid_lin, linear[m], "--", color=col, lw=1.3, alpha=0.8) + # markers show newton's discrete per-pump equilibria + ax.plot(n_cols, newton[m], ".-", color=col, lw=1.8, ms=5, label=f"mode {m}") + if xmax is not None: + ax.set_xlim(0, xmax) + top = max( + (newton[m][n_cols <= xmax].max() if np.any(n_cols <= xmax) else 0) for m in active + ) + ax.set_ylim(0, 1.15 * max(top, 1e-9)) + ax.set_xlabel("pump $D_0$") + ax.set_ylabel("modal intensity") + if active: + ax.legend(fontsize=8) + + +def main(): + graph = build_chaotic_ring() + passive = find_passive_modes(graph, method="contour") + if len(passive) == 0: + raise SystemExit("no passive modes found") + graph.graph["params"]["pump"] = np.array( + [1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()] + ) + trajectories = pump_trajectories(passive, graph, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, graph) + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + threshold_modes = tdf["threshold_lasing_modes"].to_numpy() + print("modes below the max pump (id, k, threshold, participation):") + for i in range(len(tdf)): + if thr[i] < D0_MAX: + k = from_complex(threshold_modes[i])[0] + part = _participation(threshold_modes[i], graph) + print(f" {i}: k={k:.3f} thr={thr[i]:.3f} participation={part:.1f}") + + competition = compute_mode_competition_matrix(graph, tdf) + n_modes = len(tdf) + first = float(thr[thr < np.inf].min()) + + # Full range. Sample linear on the same uniform grid newton uses: the + # event-driven sweep otherwise only stores points at mode thresholds, which + # here all cluster near 0.02 and collapse to a 2-point grid. linear is exact + # between events, so endpoint-sampling a uniform grid simply draws the line. + grid = np.linspace(first, D0_MAX, D0_STEPS) + linear = _linear_on_grid(tdf, competition, grid, n_modes) + n_cols, newton = _curves( + compute_modal_intensities_full_salt_newton(graph, tdf.copy(), D0_MAX, D0_steps=D0_STEPS) + ) + + # Zoom on the onset: the four thresholds sit in 0.016--0.026, so re-sample a + # fine grid over a small pump window (a newton run with a small max_pump is + # just a dense linspace there) to see each mode switch on in turn. + grid_z = np.linspace(first, D0_ZOOM, D0_STEPS_ZOOM) + linear_z = _linear_on_grid(tdf, competition, grid_z, n_modes) + n_cols_z, newton_z = _curves( + compute_modal_intensities_full_salt_newton( + graph, tdf.copy(), D0_ZOOM, D0_steps=D0_STEPS_ZOOM + ) + ) + + cmap = plt.get_cmap("tab10") + fig, axes = plt.subplots(1, 3, figsize=(16, 4.7)) + _draw_geometry(axes[0], graph) + _plot_li(axes[1], grid, linear, n_cols, newton, thr, cmap) + axes[1].set_title("L--I, full range (dashed = linear, solid = newton)") + _plot_li(axes[2], grid_z, linear_z, n_cols_z, newton_z, thr, cmap, xmax=D0_ZOOM) + axes[2].set_title("zoom on onset (dotted = thresholds)") + + fig.suptitle("Single ring + random chords: genuine multimode lasing", y=1.02) + fig.tight_layout() + out = HERE / "chaotic_ring_multimode.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + + def _count(arr): + return int(np.sum(arr > 1e-2 * max(arr.max(), 1e-9))) + + print(f"linear: {_count(linear[:, -1])} lasing @max") + print(f"full_salt_newton: {_count(newton[:, -1])} lasing @max") + print(f"wrote {out}") + + cut_l = 1e-2 * max(linear[:, -1].max(), 1e-9) + cut_n = 1e-2 * max(newton[:, -1].max(), 1e-9) + mode_profile_figure( + graph, + tdf, + D0_MAX, + [m for m in range(n_modes) if linear[m, -1] > cut_l], + [m for m in range(n_modes) if newton[m, -1] > cut_n], + "chaotic ring", + HERE, + a0={m: float(newton[m, -1]) for m in range(n_modes)}, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/chaotic_ring/run.sh b/examples/chaotic_ring/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/chaotic_ring/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/chord_sweep/README.md b/examples/chord_sweep/README.md new file mode 100644 index 00000000..0898ea6e --- /dev/null +++ b/examples/chord_sweep/README.md @@ -0,0 +1,22 @@ +# Chord-count sweep — density alone does not buy lasing modes + +One 14-node ring with two leads and a growing, nested set of random chords +(6 → 10 → 14; the 6-chord set is `../chaotic_ring`). The naive expectation — +more loops, denser spectrum, more lasing modes — is falsified by the measured +sweep: + +| chords | modes in window | mean participation | linear lases | newton lases | +|---|---|---|---|---| +| 6 | 6 | 7.9 | 3 | 4 | +| 10 | 8 | 8.4 | 1 | 1 | +| 14 | 7 | 8.9 | 1 | 1 | + +Adding chords delocalises the modes (participation grows) and reshuffles the +cluster under the gain, so overlap rises and the first lasing mode clamps the +gain for the rest — winner-take-all. Multimode lasing needs *localised* modes: +compare `../chaotic_ring` (a localised 4-mode cluster), `../ring_chain` +(localisation by detuning, one mode per ring) and `../mini_buffon` (extended +disorder modes, 2 of 10 lase). + +`bash run.sh` writes `chord_sweep.png` here (takes a few minutes). Figures are +not committed; re-run to reproduce. diff --git a/examples/chord_sweep/run.py b/examples/chord_sweep/run.py new file mode 100644 index 00000000..f80ab91d --- /dev/null +++ b/examples/chord_sweep/run.py @@ -0,0 +1,182 @@ +"""Chord-count sweep: density alone does not buy lasing modes. + +The chaotic-ring family made systematic: one 14-node ring with two leads and a +growing number of hard-coded random chords (6 -> 10 -> 14, nested sets so the +sweep changes one thing at a time). The naive expectation -- more chords, more +loops, denser spectrum, more lasing modes -- is **falsified** by the measured +sweep: 6 chords lase 3 (linear) / 4 (newton) modes, while 10 and 14 chords +collapse to **single-mode** lasing despite holding more modes under the gain. + +The mechanism is in the participation column: adding chords *delocalises* +the modes (the mean participation ratio grows monotonically, 7.9 -> 8.4 -> +8.9) and reshuffles the cluster under the gain, so the spatial overlap rises +and the first lasing mode clamps the gain for the others -- winner-take-all. +Multimode lasing needs *localised* modes (weak overlap), which raw loop +density does not provide: compare +``../chaotic_ring`` (6 chords, localised cluster, 4 modes), +``../ring_chain`` (localisation by detuning, one mode per ring) and +``../mini_buffon`` (extended disorder modes, 2 of 10 lase). + +Run from this directory (writes ``chord_sweep.png``; takes tens of minutes):: + + OMP_NUM_THREADS=1 python run.py +""" + +from __future__ import annotations + +import sys +import time +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import curves, linear_on_grid, threshold_modes + +import netsalt +from netsalt.modes import ( + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + mode_on_nodes, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import from_complex + +HERE = Path(__file__).resolve().parent +M = 14 # ring nodes +# Nested chord sets: the first 6 are ../chaotic_ring's, extended twice. +CHORDS = [ + (0, 7), + (1, 5), + (2, 12), + (3, 5), + (3, 10), + (7, 11), # 6 (chaotic_ring) + (4, 9), + (6, 13), + (2, 8), + (0, 10), # -> 10 + (5, 12), + (1, 9), + (8, 13), + (4, 11), # -> 14 +] +SWEEP = [6, 10, 14] +D0_MAX = 0.5 +D0_STEPS = 12 + +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 2.81, + "gamma_perp": 0.10, + "k_min": 2.55, + "k_max": 3.25, + "alpha_min": -0.05, + "alpha_max": 0.25, + "n_workers": 1, + "n_modes_max": 60, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": D0_MAX, + "D0_steps": 14, + "dielectric_params": {"method": "uniform", "inner_value": 9.0, "outer_value": 1.0, "loss": 0.0}, +} + + +def build(n_chords): + g = nx.cycle_graph(M) + g.add_edges_from(CHORDS[:n_chords]) + g.add_edge(0, M) # lead + g.add_edge(M // 2, M + 1) # lead + pos = {i: [np.cos(2 * np.pi * i / M), np.sin(2 * np.pi * i / M)] for i in range(M)} + pos[M] = [1.6, 0.0] + pos[M + 1] = [-1.6, 0.0] + positions = np.array([pos[i] for i in range(len(g))]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 12.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def _participation(mode, graph): + """Effective number of nodes the mode lives on (1 / sum p_i^2). + + Threshold modes sit at the lasing point where the operator is exactly + singular and ARPACK shift-invert cannot factorise it; a tiny imaginary + offset lifts the singularity (same guard as the field-intensity helper). + """ + k, alpha = from_complex(mode) + nudged = [float(k), float(alpha) if abs(alpha) > 1e-7 else 1e-7] + weight = np.abs(mode_on_nodes(nudged, graph, check_quality=False)) ** 2 + prob = weight / weight.sum() + return 1.0 / np.sum(prob**2) + + +def main(): + rows = [] + print( + f"{'chords':>6} {'modes':>6} {'min dk':>8} {'':>8}" + f" {'lin lases':>9} {'nwt lases':>9} {'t_nwt':>7}" + ) + for n_chords in SWEEP: + graph = build(n_chords) + tdf = threshold_modes(graph, method="contour") + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + tms = tdf["threshold_lasing_modes"].to_numpy() + n_modes = len(tdf) + ks = np.sort([from_complex(tms[i])[0] for i in range(n_modes)]) + min_dk = float(np.diff(ks).min()) if len(ks) > 1 else np.nan + partic = float(np.mean([_participation(tms[i], graph) for i in range(n_modes)])) + + competition = compute_mode_competition_matrix(graph, tdf) + first = float(thr[thr < np.inf].min()) + grid = np.linspace(first, D0_MAX, D0_STEPS) + linear = linear_on_grid(tdf, competition, grid, n_modes) + t0 = time.perf_counter() + _, newton = curves( + compute_modal_intensities_full_salt_newton(graph, tdf.copy(), D0_MAX, D0_steps=D0_STEPS) + ) + t_newton = time.perf_counter() - t0 + peak = max(linear.max(), newton.max(), 1e-9) + n_lin = int(np.sum(linear[:, -1] > 1e-2 * peak)) + n_nwt = int(np.sum(newton[:, -1] > 1e-2 * peak)) + rows.append((n_chords, n_modes, min_dk, partic, n_lin, n_nwt, t_newton)) + print( + f"{n_chords:>6} {n_modes:>6} {min_dk:>8.4f} {partic:>8.1f}" + f" {n_lin:>9} {n_nwt:>9} {t_newton:>6.0f}s" + ) + + rows = np.array(rows) + fig, axes = plt.subplots(1, 2, figsize=(10, 4)) + axes[0].plot(rows[:, 0], rows[:, 1], "s--", color="0.5", label="modes in window") + axes[0].plot(rows[:, 0], rows[:, 4], "o-", color="tab:blue", label="linear lases") + axes[0].plot(rows[:, 0], rows[:, 5], "x-", color="tab:red", label="newton lases") + axes[0].set_xlabel("number of chords") + axes[0].set_ylabel("count") + axes[0].legend(fontsize=8) + axes[0].set_title("denser spectrum, fewer lasing modes") + axes[1].plot(rows[:, 0], rows[:, 3], "o-", color="tab:purple") + axes[1].set_xlabel("number of chords") + axes[1].set_ylabel("mean participation ratio") + axes[1].set_title("the mechanism: modes delocalise") + fig.suptitle("Ring + chords: density alone does not buy lasing modes", y=1.02) + fig.tight_layout() + out = HERE / "chord_sweep.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") + + +if __name__ == "__main__": + main() diff --git a/examples/chord_sweep/run.sh b/examples/chord_sweep/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/chord_sweep/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/dense_ring/README.md b/examples/dense_ring/README.md new file mode 100644 index 00000000..c99a6c13 --- /dev/null +++ b/examples/dense_ring/README.md @@ -0,0 +1,13 @@ +# Denser chord ring — newton-vs-linear consistency + +A 16-node ring with 10 hard-coded chords: a bigger, more strongly-competing +companion to `../chaotic_ring`. Both solvers share the onset slope at +threshold, so they coincide just above it; above threshold they diverge because +`full_salt_newton` re-solves each lasing mode's profile at the operating pump +(the spatial holes, and hence the competition, shift) while `linear` freezes +the profiles. The signature: the dominant mode tracks linear closely, while the +secondary modes are reshuffled — one lights up earlier and stronger, another +later and largely suppressed. + +`bash run.sh` writes `dense_ring_compare.png` here (and prints a per-mode +onset/slope/intensity table). Figures are not committed; re-run to reproduce. diff --git a/examples/dense_ring/run.py b/examples/dense_ring/run.py new file mode 100644 index 00000000..cd597c9f --- /dev/null +++ b/examples/dense_ring/run.py @@ -0,0 +1,221 @@ +"""Newton vs linear on a *denser* chord ring: consistency of the L--I curves. + +A companion to ``../chaotic_ring/run.py`` that checks the operator-level +``full_salt_newton`` solver against the near-threshold ``linear`` model on a +bigger, more strongly-competing graph: a **16-node ring with 10 random chords** +(vs 14 nodes / 6 chords). The question this answers is *what should differ between +the two solvers, and why*: + +* Both are built to share the same onset slope ``1/(T_μμ·D0_thr)`` at each mode's + threshold, so just above threshold they **coincide**. +* ``linear`` then freezes every mode's profile at its own threshold and uses a + fixed competition matrix ``T`` -> strictly piecewise-linear curves. +* ``full_salt_newton`` re-solves each lasing mode's ``(k_μ, a_μ)`` *and spatial + profile* at the operating pump, under the shared saturated (hole-burnt) + operator. As the pump rises the holes deepen and the profiles/overlaps shift, so + the effective competition changes. + +The visible consequences, printed and plotted below: + +* the **dominant** mode tracks the linear curve closely (same onset, near-equal + initial slope), then its slope drifts as saturation sets in; +* the **secondary** modes are *reshuffled* -- they switch on at different pumps and + reach different intensities than the frozen-profile linear model predicts (here + one secondary lights up earlier and stronger, another later and largely + suppressed). That reshuffling is the operator-level mode interaction the linear + model cannot see. + +The chord layout is hard-coded (from a small seed scan) so the result is +reproducible regardless of the NumPy RNG. + +Run:: + + OMP_NUM_THREADS=1 python run.py +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import mode_profile_figure + +import netsalt +from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + find_passive_modes, + find_threshold_lasing_modes, + pump_trajectories, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import from_complex + +HERE = Path(__file__).resolve().parent +M = 16 # ring nodes +# 10 chords across the ring + a lead on node 0 and node M // 2 (hard-coded so the +# spectrum is reproducible regardless of the NumPy RNG). +CHORDS = [(0, 13), (1, 14), (2, 10), (3, 8), (4, 9), (6, 13), (7, 15), (9, 12), (10, 15), (11, 15)] + +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 2.9, + "gamma_perp": 0.22, + "k_min": 2.55, + "k_max": 3.25, + "alpha_min": -0.05, + "alpha_max": 0.25, + "n_workers": 1, + "n_modes_max": 50, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 0.5, + "D0_steps": 14, + "dielectric_params": {"method": "uniform", "inner_value": 9.0, "outer_value": 1.0, "loss": 0.0}, +} +D0_MAX = 0.5 +D0_STEPS = 24 + + +def build_dense_ring(): + """A 16-node ring with 10 chords and two leads (open cavity).""" + g = nx.cycle_graph(M) + g.add_edges_from(CHORDS) + g.add_edge(0, M) # lead + g.add_edge(M // 2, M + 1) # lead + pos = {i: [np.cos(2 * np.pi * i / M), np.sin(2 * np.pi * i / M)] for i in range(M)} + pos[M] = [1.7, 0.0] + pos[M + 1] = [-1.7, 0.0] + positions = np.array([pos[i] for i in range(len(g))]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 12.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def _endpoint(df): + cols = sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + return np.nan_to_num(df[("modal_intensities", cols[-1])].to_numpy(dtype=float)) + + +def _curves(df): + cols = np.array( + sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + ) + return cols, np.nan_to_num(df[[("modal_intensities", c) for c in cols]].to_numpy(dtype=float)) + + +def _onset_slope(pumps, y): + """(onset pump, near-threshold slope) -- slope fit just above the turn-on.""" + idx = np.where(y > 1e-3)[0] + if idx.size < 3: + return np.nan, np.nan + seg = idx[1 : min(idx.size, 5)] # skip the activation point itself + slope = np.polyfit(pumps[seg], y[seg], 1)[0] + return float(pumps[idx[0]]), float(slope) + + +def _draw_geometry(ax, graph): + pos = {n: graph.nodes[n]["position"] for n in graph.nodes} + chord_set = {tuple(sorted(c)) for c in CHORDS} + for u, v in graph.edges(): + x = [pos[u][0], pos[v][0]] + y = [pos[u][1], pos[v][1]] + su, sv = min(u, v), max(u, v) + if sv >= M: + ax.plot(x, y, color="crimson", lw=2.0, ls="--", zorder=1) + elif (su, sv) in chord_set: + ax.plot(x, y, color="royalblue", lw=2.0, zorder=2) + else: + ax.plot(x, y, color="0.6", lw=2.0, zorder=1) + for n in graph.nodes: + ax.scatter(*pos[n], s=130 if n < M else 100, color="white", edgecolor="black", zorder=3) + ax.set_aspect("equal") + ax.axis("off") + ax.set_title(f"{M}-node ring + {len(CHORDS)} chords") + + +def main(): + graph = build_dense_ring() + passive = find_passive_modes(graph, method="contour") + if len(passive) == 0: + raise SystemExit("no passive modes found") + graph.graph["params"]["pump"] = np.array( + [1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()] + ) + trajectories = pump_trajectories(passive, graph, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, graph) + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + threshold_modes = tdf["threshold_lasing_modes"].to_numpy() + n_modes = len(tdf) + + competition = compute_mode_competition_matrix(graph, tdf) + first = float(thr[thr < np.inf].min()) + grid = np.linspace(first, D0_MAX, D0_STEPS) + linear = np.zeros((n_modes, grid.size)) + for j, d0 in enumerate(grid): + linear[:, j] = _endpoint(compute_modal_intensities(tdf.copy(), d0, competition)) + n_cols, newton = _curves( + compute_modal_intensities_full_salt_newton(graph, tdf.copy(), D0_MAX, D0_steps=D0_STEPS) + ) + + peak = max(linear.max(), newton.max(), 1e-9) + active = [m for m in range(n_modes) if max(linear[m].max(), newton[m].max()) > 1e-2 * peak] + print("mode k thr onset(lin/nwt) slope(lin/nwt) I@max(lin/nwt)") + for m in active: + k = from_complex(threshold_modes[m])[0] + on_l, sl_l = _onset_slope(grid, linear[m]) + on_n, sl_n = _onset_slope(n_cols, newton[m]) + print( + f"{m:>3} {k:.3f} {thr[m]:.3f} {on_l:.3f} / {on_n:.3f} " + f"{sl_l:7.1f} / {sl_n:7.1f} {linear[m, -1]:7.2f} / {newton[m, -1]:7.2f}" + ) + + cmap = plt.get_cmap("tab10") + fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + _draw_geometry(axes[0], graph) + for m in active: + col = cmap(m % 10) + axes[1].plot(grid, linear[m], "--", color=col, lw=1.3, alpha=0.8) + axes[1].plot(n_cols, newton[m], ".-", color=col, lw=1.8, ms=4, label=f"mode {m}") + axes[1].set_xlabel("pump $D_0$") + axes[1].set_ylabel("modal intensity") + axes[1].set_title("dashed = linear, solid = newton") + axes[1].legend(fontsize=8) + fig.suptitle("Denser chord ring: full_salt_newton vs linear", y=1.0) + fig.tight_layout() + cut = 1e-2 * peak + mode_profile_figure( + graph, + tdf, + D0_MAX, + [m for m in range(n_modes) if linear[m, -1] > cut], + [m for m in range(n_modes) if newton[m, -1] > cut], + "dense chord ring", + HERE, + a0={m: float(newton[m, -1]) for m in range(n_modes)}, + ) + + out = HERE / "dense_ring_compare.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") + + +if __name__ == "__main__": + main() diff --git a/examples/dense_ring/run.sh b/examples/dense_ring/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/dense_ring/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/line_PRA/README.md b/examples/line_PRA/README.md new file mode 100644 index 00000000..86678fe4 --- /dev/null +++ b/examples/line_PRA/README.md @@ -0,0 +1,60 @@ +# line_PRA — the Ge–Chong–Stone 1D slab laser + +This example is the partially-pumped 1D slab resonator of Ge, Chong & Stone, +*Steady-state ab initio laser theory: generalizations and analytic results*, +[PRA **82**, 063824 (2010)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.82.063824) +([arXiv:1008.0628](https://arxiv.org/abs/1008.0628)), Figs. 3/5/6: + +- cavity of length 1, open on both sides; +- refractive index n = 1.5 on 0 < x < 0.25 and n = 3 on 0.25 < x < 1; +- gain centre k_a = 15, width gamma_perp = 3; +- only the left half (0 < x < 0.5) pumped. + +`create_graph.py` builds it as a 10-edge line graph (one n = 1 lead edge at +each end carrying the open boundary, 8 inner edges of length 0.125). + +## Run the example + +```bash +bash run.sh # create the graph, then python -m netsalt lasing config.yaml +``` + +## Reproduce the Fig. 6 validation + +The paper's Fig. 6 shows modal intensity vs pump for two lasing modes, both as +the exact SALT (numerical solution of their Eq. 28, open symbols) and as the +single-pole approximation (SPA, solid lines). netsalt's `linear` solver *is* +the SPA, and `full_salt_newton` is an operator-level SALT that should track +the exact symbols (dominant mode suppressed below the SPA once the second mode +turns on, second mode above it). + +```bash +python create_graph.py # graph.json / index.yaml / pump.yaml +OMP_NUM_THREADS=1 python compare_to_pra_fig6.py +``` + +This runs the pipeline (cached in `out/`, so re-runs are fast), overlays both +solvers on the digitized Fig. 6 data, writes +`figures/pra_fig6_compare.pdf`, and prints a checkpoint table. Expected +agreement (D0 = 1.258, the figure edge): `linear` vs paper SPA dominant to +< 1 %; `newton` vs paper exact dominant to ~1 %; second mode within ~10 %, +dominated by a +0.03 offset of its interacting threshold (paper exact 0.892, +SPA 0.899; netsalt ~0.92–0.93). That offset traces to the second mode's +*noninteracting* threshold being ~0.3 % above the paper's Fig. 3(b) value +(0.6641 vs 0.662) — netsalt matches the four modes at and below the gain +centre to < 0.1 % but sits +0.3–0.4 % high on the two above it — amplified +~6x by gain-clamping proximity in the interacting-threshold formula. It is +converged in netsalt (identical at 2–4x finer pump stepping), so it is a +sub-half-percent model difference (netsalt's quantum-graph dispersion vs the +paper's 20-state TCF basis), not a numerics artifact of either solver. + +`data/ge_fig6_digitized.csv` is the digitized reference data. Regenerate it +from the paper itself with + +```bash +python digitize_pra_fig6.py # downloads the arXiv PDF; needs poppler + scipy +``` + +(axis calibration from the tick labels is good to ~0.3 %; in the single-mode +regime the exact symbols ride on the SPA lines, so only the line is recovered +there and the `spa_*` series carry small bumps where markers merged into it). diff --git a/examples/line_PRA/compare_to_pra_fig6.py b/examples/line_PRA/compare_to_pra_fig6.py new file mode 100644 index 00000000..ff6bda21 --- /dev/null +++ b/examples/line_PRA/compare_to_pra_fig6.py @@ -0,0 +1,159 @@ +"""Compare netsalt's modal-intensity solvers to Ge-Chong-Stone Fig. 6. + +The ``line_PRA`` graph is exactly the 1D slab of Ge-Chong-Stone, PRA 82, +063824 (arXiv:1008.0628), Figs. 3/5/6: length 1, n = 1.5 on x < 0.25 and n = 3 +on x > 0.25, open on both sides, gain k_a = 15, gamma_perp = 3, left half +pumped. This script runs the shared netsalt pipeline (cached in ``out/``), then +overlays + +* ``linear`` -- netsalt's competition matrix; should coincide with the + paper's single-pole approximation (SPA) lines; +* ``full_salt_newton`` -- the operator-level SALT; should track the paper's + exact symbols (Eq. 28): dominant suppressed below the + SPA, second mode above it, once the second mode is on; + +against ``data/ge_fig6_digitized.csv`` (regenerate with +``digitize_pra_fig6.py``). Writes ``figures/pra_fig6_compare.pdf`` and prints a +checkpoint table. Run from this directory (after ``python create_graph.py`` or +``bash run.sh``):: + + OMP_NUM_THREADS=1 python compare_to_pra_fig6.py +""" + +import csv +import os +from collections import defaultdict +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +HERE = Path(__file__).resolve().parent +D0_MAX = 1.27 # the Fig. 6 pump range + + +def load_reference(): + series = defaultdict(list) + with open(HERE / "data" / "ge_fig6_digitized.csv") as f: + for row in csv.DictReader(r for r in f if not r.startswith("#")): + series[row["series"]].append((float(row["D0"]), float(row["intensity"]))) + return {k: np.array(v) for k, v in series.items()} + + +def run_netsalt(): + from netsalt.config_loader import load_config + from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + ) + from netsalt.pipeline import ( + _attach_pump_to_graph, + step_compute_mode_trajectories, + step_create_pump_profile, + step_create_quantum_graph, + step_find_passive_modes, + step_find_threshold_modes, + step_scan_frequencies, + ) + + p = load_config(HERE / "config.yaml") + ids = p.get("lasing_modes_id") + qg = step_create_quantum_graph(p) + qualities = step_scan_frequencies(p, qg) + passive = step_find_passive_modes(p, qg, qualities) + pump = step_create_pump_profile(p, qg, passive, ids) + traj = step_compute_mode_trajectories(p, qg, passive, pump, ids) + tdf = step_find_threshold_modes(p, qg, traj, pump, ids) + qg = _attach_pump_to_graph(p, qg, pump) + + T = compute_mode_competition_matrix(qg, tdf) + lin = compute_modal_intensities(tdf.copy(), D0_MAX, T) + new = compute_modal_intensities_full_salt_newton(qg, tdf.copy(), D0_MAX, D0_steps=16) + return tdf, lin, new + + +def curves(df): + cols = [c for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + pumps = np.array(sorted(c[1] for c in cols), dtype=float) + data = np.nan_to_num(df[[("modal_intensities", c) for c in pumps]].to_numpy(float)) + return pumps, data + + +def main(): + os.chdir(HERE) + ref = load_reference() + tdf, lin, new = run_netsalt() + + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + pl, dl = curves(lin) + pn, dn = curves(new) + lasing = [i for i in np.argsort(thresholds) if dn[i].max() > 0][:2] + dom, sec = lasing[0], (lasing[1] if len(lasing) > 1 else None) + + # one color per mode; series are told apart by style (thin solid = paper SPA, + # open markers = paper exact, dashed = netsalt linear, x = netsalt newton) + c_dom, c_sec = "firebrick", "steelblue" + fig, ax = plt.subplots(figsize=(6, 4.5)) + ax.plot(*ref["spa_dominant"].T, "-", c=c_dom, lw=1, alpha=0.6, label="paper SPA") + # below its threshold the digitized "line" is the chain of zero-intensity + # exact markers; show the SPA line only where the mode is on + spa2 = ref["spa_second"][ref["spa_second"][:, 1] > 0.008] + slope, icpt = np.polyfit(spa2[:, 0], spa2[:, 1], 1) + spa2 = np.vstack([[-icpt / slope, 0.0], spa2]) # extend down to its zero crossing + ax.plot(*spa2.T, "-", c=c_sec, lw=1, alpha=0.6) + ax.plot(*ref["exact_dominant"].T, "s", mfc="none", c=c_dom, ms=6, label="paper exact (Eq. 28)") + ax.plot(*ref["exact_second"].T, "o", mfc="none", c=c_sec, ms=6) + ax.plot(pl, dl[dom], "--", c=c_dom, lw=2, label=f"netsalt linear (mode {dom})") + ax.plot(pn, dn[dom], "x-", c=c_dom, lw=1, ms=7, label=f"netsalt newton (mode {dom})") + if sec is not None: + ax.plot(pl, dl[sec], "--", c=c_sec, lw=2, label=f"netsalt linear (mode {sec})") + ax.plot(pn, dn[sec], "x-", c=c_sec, lw=1, ms=7, label=f"netsalt newton (mode {sec})") + ax.set_xlabel(r"pump strength $D_0$") + ax.set_ylabel(r"modal intensity $I_\mu$") + ax.set_xlim(0.58, 1.3) + ax.set_ylim(0, 0.33) + ax.legend(loc="upper left", fontsize=8) + ax.set_title("line_PRA vs Ge-Chong-Stone PRA 82, 063824 Fig. 6") + fig.tight_layout() + out = HERE / "figures" / "pra_fig6_compare.pdf" + out.parent.mkdir(exist_ok=True) + fig.savefig(out) + print(f"wrote {out}") + + # checkpoint table at the figure edge + d0c = 1.258 + + def at(p, row): + on = row > 0 + return float(np.interp(d0c, p[on], row[on])) if on.sum() > 1 else 0.0 + + def ref_at(name): + d = ref[name] + return float(np.interp(d0c, d[:, 0], d[:, 1])) + + print(f"\ncheckpoint at D0 = {d0c}:") + print( + f" dominant: paper exact {ref_at('exact_dominant'):.3f} newton {at(pn, dn[dom]):.3f}" + f" | paper SPA {ref_at('spa_dominant'):.3f} linear {at(pl, dl[dom]):.3f}" + ) + if sec is not None: + print( + f" second : paper exact {ref_at('exact_second'):.3f} newton {at(pn, dn[sec]):.3f}" + f" | paper SPA {ref_at('spa_second'):.3f} linear {at(pl, dl[sec]):.3f}" + ) + print(f" first threshold: netsalt {thresholds[dom]:.4f} (paper ~0.61)") + if sec is not None: + ithr = np.asarray(new["interacting_lasing_thresholds"]).ravel() + print( + f" 2nd interacting threshold: linear " + f"{np.asarray(lin['interacting_lasing_thresholds']).ravel()[sec]:.4f}, " + f"newton <= {ithr[sec]:.4f} (paper exact 0.892, SPA 0.899)" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/line_PRA/config.yaml b/examples/line_PRA/config.yaml index 4047b13c..572a5596 100644 --- a/examples/line_PRA/config.yaml +++ b/examples/line_PRA/config.yaml @@ -27,6 +27,11 @@ D0_max: 2.0 D0_steps: 30 intensities_D0_max: 4.0 +# Modal-intensity solver (issue #42): "linear" (default, fast near-threshold +# SALT) or "full_salt_newton" (operator-level SALT, precise above threshold). +# See doc/source/lasing.rst and benchmark/bench_salt.py. +intensity_method: linear + pump_mode: custom pump_custom_path: pump.yaml diff --git a/examples/line_PRA/data/.gitignore b/examples/line_PRA/data/.gitignore new file mode 100644 index 00000000..d6cabc9a --- /dev/null +++ b/examples/line_PRA/data/.gitignore @@ -0,0 +1,2 @@ +arxiv_1008.0628.pdf +ge_fig6-*.png diff --git a/examples/line_PRA/data/ge_fig6_digitized.csv b/examples/line_PRA/data/ge_fig6_digitized.csv new file mode 100644 index 00000000..4e7d0f74 --- /dev/null +++ b/examples/line_PRA/data/ge_fig6_digitized.csv @@ -0,0 +1,273 @@ +# Ge-Chong-Stone PRA 82, 063824 (arXiv:1008.0628) Fig. 6, digitized +# exact_* : symbols (numerical solution of Eq. 28, exact SALT) +# spa_* : solid lines (single-pole approximation, Eqs. 40, 44-45) +# NB: in the single-mode regime the symbols ride on the lines, so the +# spa_* series carry small bumps where merged markers were not separable. +series,D0,intensity +exact_dominant,1.0640,0.1598 +exact_dominant,1.0887,0.1661 +exact_dominant,1.1133,0.1722 +exact_dominant,1.1384,0.1786 +exact_dominant,1.1629,0.1850 +exact_dominant,1.1876,0.1913 +exact_dominant,1.2125,0.1975 +exact_dominant,1.2372,0.2041 +exact_dominant,1.2582,0.2104 +exact_second,1.0392,0.0429 +exact_second,1.0641,0.0501 +exact_second,1.0887,0.0574 +exact_second,1.1134,0.0647 +exact_second,1.1382,0.0721 +exact_second,1.1629,0.0796 +exact_second,1.1876,0.0870 +exact_second,1.2123,0.0945 +exact_second,1.2370,0.1020 +exact_second,1.2579,0.1096 +exact_total,1.2584,0.3200 +exact_total,0.9010,0.1238 +exact_total,0.9053,0.1251 +exact_total,0.9096,0.1283 +exact_total,0.9138,0.1326 +exact_total,0.9181,0.1328 +exact_total,0.9224,0.1340 +exact_total,0.9266,0.1363 +exact_total,0.9309,0.1385 +exact_total,0.9352,0.1415 +exact_total,0.9395,0.1451 +exact_total,0.9437,0.1443 +exact_total,0.9480,0.1413 +exact_total,0.9523,0.1496 +exact_total,0.9565,0.1525 +exact_total,0.9608,0.1553 +exact_total,0.9651,0.1590 +exact_total,0.9693,0.1568 +exact_total,0.9736,0.1543 +exact_total,0.9779,0.1630 +exact_total,0.9822,0.1663 +exact_total,0.9864,0.1693 +exact_total,0.9907,0.1718 +exact_total,0.9950,0.1696 +exact_total,0.9992,0.1733 +exact_total,1.0035,0.1763 +exact_total,1.0078,0.1803 +exact_total,1.0120,0.1831 +exact_total,1.0163,0.1846 +exact_total,1.0206,0.1823 +exact_total,1.0248,0.1873 +exact_total,1.0291,0.1903 +exact_total,1.0334,0.1944 +exact_total,1.0377,0.1971 +exact_total,1.0419,0.1973 +exact_total,1.0462,0.1951 +exact_total,1.0505,0.2009 +exact_total,1.0547,0.2041 +exact_total,1.0590,0.2088 +exact_total,1.0633,0.2136 +exact_total,1.0675,0.2103 +exact_total,1.0718,0.2081 +exact_total,1.0761,0.2143 +exact_total,1.0804,0.2188 +exact_total,1.0846,0.2233 +exact_total,1.0889,0.2249 +exact_total,1.0932,0.2233 +exact_total,1.0974,0.2226 +exact_total,1.1017,0.2276 +exact_total,1.1060,0.2331 +exact_total,1.1102,0.2378 +exact_total,1.1145,0.2404 +exact_total,1.1188,0.2358 +exact_total,1.1231,0.2379 +exact_total,1.1273,0.2409 +exact_total,1.1316,0.2473 +exact_total,1.1359,0.2523 +exact_total,1.1401,0.2533 +exact_total,1.1444,0.2486 +exact_total,1.1487,0.2519 +exact_total,1.1529,0.2551 +exact_total,1.1572,0.2611 +exact_total,1.1615,0.2673 +exact_total,1.1658,0.2659 +exact_total,1.1700,0.2619 +exact_total,1.1743,0.2656 +exact_total,1.1786,0.2689 +exact_total,1.1828,0.2759 +exact_total,1.1871,0.2749 +exact_total,1.1914,0.2787 +exact_total,1.1956,0.2756 +exact_total,1.1999,0.2789 +exact_total,1.2042,0.2874 +exact_total,1.2085,0.2887 +exact_total,1.2127,0.2887 +exact_total,1.2170,0.2916 +exact_total,1.2213,0.2892 +exact_total,1.2255,0.2922 +exact_total,1.2298,0.3021 +exact_total,1.2341,0.3027 +exact_total,1.2383,0.3029 +exact_total,1.2426,0.3041 +exact_total,1.2469,0.3029 +exact_total,1.2512,0.3056 +exact_total,1.2554,0.3079 +exact_total,1.2597,0.3101 +spa_dominant,0.6121,0.0039 +spa_dominant,0.6183,0.0050 +spa_dominant,0.6244,0.0040 +spa_dominant,0.6306,0.0065 +spa_dominant,0.6396,0.0102 +spa_dominant,0.6458,0.0128 +spa_dominant,0.6548,0.0170 +spa_dominant,0.6610,0.0197 +spa_dominant,0.6700,0.0235 +spa_dominant,0.6795,0.0277 +spa_dominant,0.6856,0.0303 +spa_dominant,0.6946,0.0340 +spa_dominant,0.7041,0.0382 +spa_dominant,0.7103,0.0408 +spa_dominant,0.7193,0.0447 +spa_dominant,0.7288,0.0488 +spa_dominant,0.7350,0.0515 +spa_dominant,0.7445,0.0555 +spa_dominant,0.7535,0.0593 +spa_dominant,0.7596,0.0620 +spa_dominant,0.7691,0.0662 +spa_dominant,0.7781,0.0700 +spa_dominant,0.7843,0.0727 +spa_dominant,0.7933,0.0767 +spa_dominant,0.8023,0.0805 +spa_dominant,0.8085,0.0830 +spa_dominant,0.8175,0.0870 +spa_dominant,0.8265,0.0908 +spa_dominant,0.8327,0.0935 +spa_dominant,0.8417,0.0973 +spa_dominant,0.8507,0.1013 +spa_dominant,0.8569,0.1040 +spa_dominant,0.8659,0.1078 +spa_dominant,0.8749,0.1116 +spa_dominant,0.8811,0.1143 +spa_dominant,0.8901,0.1183 +spa_dominant,0.8996,0.1223 +spa_dominant,0.9072,0.1241 +spa_dominant,0.9281,0.1305 +spa_dominant,0.9489,0.1363 +spa_dominant,0.9551,0.1381 +spa_dominant,0.9760,0.1440 +spa_dominant,0.9822,0.1456 +spa_dominant,1.0026,0.1513 +spa_dominant,1.0230,0.1571 +spa_dominant,1.0291,0.1588 +spa_dominant,1.0495,0.1645 +spa_dominant,1.0557,0.1663 +spa_dominant,1.0619,0.1680 +spa_dominant,1.0680,0.1696 +spa_dominant,1.0742,0.1715 +spa_dominant,1.0804,0.1731 +spa_dominant,1.0865,0.1748 +spa_dominant,1.0927,0.1766 +spa_dominant,1.0989,0.1783 +spa_dominant,1.1050,0.1800 +spa_dominant,1.1112,0.1818 +spa_dominant,1.1174,0.1835 +spa_dominant,1.1235,0.1851 +spa_dominant,1.1297,0.1870 +spa_dominant,1.1359,0.1886 +spa_dominant,1.1420,0.1903 +spa_dominant,1.1482,0.1921 +spa_dominant,1.1544,0.1938 +spa_dominant,1.1605,0.1956 +spa_dominant,1.1667,0.1973 +spa_dominant,1.1729,0.1989 +spa_dominant,1.1790,0.2008 +spa_dominant,1.1852,0.2024 +spa_dominant,1.1914,0.2041 +spa_dominant,1.1975,0.2059 +spa_dominant,1.2037,0.2076 +spa_dominant,1.2099,0.2093 +spa_dominant,1.2160,0.2111 +spa_dominant,1.2222,0.2128 +spa_dominant,1.2284,0.2144 +spa_dominant,1.2345,0.2163 +spa_dominant,1.2407,0.2179 +spa_dominant,1.2469,0.2196 +spa_dominant,1.2531,0.2214 +spa_dominant,1.2592,0.2231 +spa_second,0.6168,0.0014 +spa_second,0.6235,0.0014 +spa_second,0.6301,0.0014 +spa_second,0.6368,0.0027 +spa_second,0.6448,0.0053 +spa_second,0.6529,0.0014 +spa_second,0.6595,0.0014 +spa_second,0.6676,0.0055 +spa_second,0.6757,0.0027 +spa_second,0.6823,0.0014 +spa_second,0.6899,0.0048 +spa_second,0.6965,0.0048 +spa_second,0.7046,0.0014 +spa_second,0.7112,0.0029 +spa_second,0.7188,0.0053 +spa_second,0.7264,0.0015 +spa_second,0.7331,0.0014 +spa_second,0.7411,0.0053 +spa_second,0.7478,0.0039 +spa_second,0.7554,0.0014 +spa_second,0.7634,0.0047 +spa_second,0.7701,0.0052 +spa_second,0.7777,0.0014 +spa_second,0.7843,0.0024 +spa_second,0.7924,0.0055 +spa_second,0.8000,0.0022 +spa_second,0.8066,0.0014 +spa_second,0.8147,0.0053 +spa_second,0.8213,0.0043 +spa_second,0.8294,0.0014 +spa_second,0.8374,0.0047 +spa_second,0.8441,0.0052 +spa_second,0.8522,0.0014 +spa_second,0.8588,0.0025 +spa_second,0.8669,0.0055 +spa_second,0.8749,0.0014 +spa_second,0.8816,0.0014 +spa_second,0.8892,0.0053 +spa_second,0.8958,0.0042 +spa_second,0.9039,0.0017 +spa_second,0.9276,0.0068 +spa_second,0.9518,0.0128 +spa_second,0.9760,0.0187 +spa_second,1.0002,0.0245 +spa_second,1.0244,0.0303 +spa_second,1.0310,0.0320 +spa_second,1.0377,0.0335 +spa_second,1.0443,0.0352 +spa_second,1.0509,0.0368 +spa_second,1.0576,0.0385 +spa_second,1.0642,0.0400 +spa_second,1.0709,0.0417 +spa_second,1.0775,0.0432 +spa_second,1.0842,0.0448 +spa_second,1.0908,0.0465 +spa_second,1.0974,0.0482 +spa_second,1.1041,0.0497 +spa_second,1.1107,0.0513 +spa_second,1.1174,0.0528 +spa_second,1.1240,0.0545 +spa_second,1.1306,0.0562 +spa_second,1.1373,0.0578 +spa_second,1.1439,0.0593 +spa_second,1.1506,0.0610 +spa_second,1.1572,0.0625 +spa_second,1.1639,0.0642 +spa_second,1.1705,0.0658 +spa_second,1.1771,0.0675 +spa_second,1.1838,0.0690 +spa_second,1.1904,0.0707 +spa_second,1.1971,0.0722 +spa_second,1.2037,0.0738 +spa_second,1.2104,0.0755 +spa_second,1.2170,0.0772 +spa_second,1.2236,0.0787 +spa_second,1.2303,0.0803 +spa_second,1.2369,0.0818 +spa_second,1.2436,0.0835 +spa_second,1.2502,0.0852 +spa_second,1.2568,0.0868 +spa_second,1.2635,0.0883 diff --git a/examples/line_PRA/digitize_pra_fig6.py b/examples/line_PRA/digitize_pra_fig6.py new file mode 100644 index 00000000..2d3a7027 --- /dev/null +++ b/examples/line_PRA/digitize_pra_fig6.py @@ -0,0 +1,165 @@ +"""Digitize Fig. 6 of Ge-Chong-Stone, PRA 82, 063824 (arXiv:1008.0628). + +Regenerates ``data/ge_fig6_digitized.csv``, the reference data that +``compare_to_pra_fig6.py`` overlays on the netsalt solvers. Open symbols in the +figure are the exact SALT (numerical solution of the paper's Eq. 28), solid +lines the single-pole approximation (Eqs. 40, 44-45): red squares = dominant +mode, blue circles = second mode, black triangles = total. + +Method: download the arXiv PDF (if not already next to this script), render the +Fig. 6 region at 600 dpi with ``pdftoppm`` (poppler), color-separate the +markers from the lines, and calibrate the axes from the tick-label text +centroids (MATLAB centers labels on ticks; calibration residual ~0.3%). +Isolated blobs of a color are markers; the long thin component is the SPA line +-- markers riding *on* the line (single-mode regime) merge into it, so there +the ``spa_*`` series carries small bumps and the exact symbols are only +recovered where they separate from the line. The crop box is pinned to the v1 +arXiv PDF (Fig. 6 is on page 12). + +Requires: poppler (``pdftoppm``, ``pdftotext``), scipy, matplotlib. +""" + +import subprocess +import urllib.request +from pathlib import Path + +import numpy as np +from matplotlib.image import imread +from scipy import ndimage + +HERE = Path(__file__).resolve().parent +PDF = HERE / "data" / "arxiv_1008.0628.pdf" +OUT = HERE / "data" / "ge_fig6_digitized.csv" +ARXIV_URL = "https://arxiv.org/pdf/1008.0628" +# Fig. 6 crop on its page, in 600-dpi pixels (pinned to the v1 arXiv PDF) +CROP = {"x": 2900, "y": 300, "W": 2200, "H": 1500} + + +def fig6_page(pdf): + """Page number holding the Fig. 6 caption.""" + txt = subprocess.run(["pdftotext", str(pdf), "-"], capture_output=True, text=True).stdout + for i, page in enumerate(txt.split("\f"), 1): + if "FIG. 6" in page: + return i + raise RuntimeError("FIG. 6 not found in PDF text") + + +def render(pdf, page, png_prefix): + subprocess.run( + ["pdftoppm", "-png", "-r", "600", "-f", str(page), "-l", str(page)] + + ["-x", str(CROP["x"]), "-y", str(CROP["y"]), "-W", str(CROP["W"]), "-H", str(CROP["H"])] + + [str(pdf), str(png_prefix)], + check=True, + ) + (out,) = png_prefix.parent.glob(png_prefix.name + "*.png") + return out + + +def blob_centroids(mask, min_pix=40): + lab, n = ndimage.label(mask) + out = [] + for i in range(1, n + 1): + ys, xs = np.where(lab == i) + if len(xs) >= min_pix: + out.append((xs.mean(), ys.mean())) + return out + + +def cluster_1d(vals, gap): + vals = sorted(vals) + groups, cur = [], [vals[0]] + for v in vals[1:]: + if v - cur[-1] <= gap: + cur.append(v) + else: + groups.append(float(np.mean(cur))) + cur = [v] + groups.append(float(np.mean(cur))) + return groups + + +def digitize(png): + img = imread(png)[..., :3] + H, _W = img.shape[:2] + r, g, b = img[..., 0], img[..., 1], img[..., 2] + red = (r > 0.6) & (g < 0.35) & (b < 0.35) + blue = (b > 0.6) & (r < 0.35) & (g < 0.35) + black = (r < 0.25) & (g < 0.25) & (b < 0.25) + + # plot frame = longest black vertical/horizontal runs + vert = np.where(black.sum(0) > 0.5 * H)[0] + x_left, x_right = int(vert.min()), int(vert.max()) + horiz = np.where(black[:, x_left:x_right].sum(1) > 0.4 * (x_right - x_left))[0] + y_top, y_bot = int(horiz.min()), int(horiz.max()) + + # axis calibration from the tick-label text centroids + # label band only (the axis title sits lower and must not join the clusters) + xs_lab = cluster_1d([c[0] for c in blob_centroids(black[y_bot + 5 : y_bot + 90, :])], gap=60) + assert len(xs_lab) == 4, f"expected x labels 0.6 0.8 1.0 1.2, got {len(xs_lab)}" + cx = np.polyfit(xs_lab, [0.6, 0.8, 1.0, 1.2], 1) + left = black[:, max(0, x_left - 130) : x_left - 8] + ys_lab = cluster_1d([c[1] for c in blob_centroids(left) if c[1] < y_bot + 40], gap=40) + assert len(ys_lab) == 4, f"expected y labels 0.3 0.2 0.1 0, got {len(ys_lab)}" + cy = np.polyfit(ys_lab, [0.3, 0.2, 0.1, 0.0], 1) + + inside = np.zeros_like(red) + inside[y_top + 3 : y_bot - 2, x_left + 3 : x_right - 2] = True + + def extract(mask): + lab, n = ndimage.label(mask & inside) + markers, line_pts = [], [] + for i in range(1, n + 1): + ys, xs = np.where(lab == i) + w, h = np.ptp(xs) + 1, np.ptp(ys) + 1 + if w > 120: # the SPA line (with any markers riding on it merged in) + # a column crossing a hollow marker outline is ~3x thicker than + # the bare line; skip those so the line series stays smooth + cols, counts = np.unique(xs, return_counts=True) + clean = cols[counts <= 1.5 * np.median(counts)] + line_pts.extend((c, np.median(ys[xs == c])) for c in clean) + elif 12 < w < 90 and 12 < h < 90: # an isolated open marker + markers.append((xs.mean(), ys.mean())) + return np.array(markers).reshape(-1, 2), np.array(line_pts).reshape(-1, 2) + + series = {} + for _color, mask, exact, spa in [ + ("red", red, "exact_dominant", "spa_dominant"), + ("blue", blue, "exact_second", "spa_second"), + ("black", black, "exact_total", "exact_total"), + ]: + markers, line = extract(mask) + for name, pts in [(exact, markers), (spa, line)]: + if not len(pts): + continue + d = np.column_stack([np.polyval(cx, pts[:, 0]), np.polyval(cy, pts[:, 1])]) + d = d[np.argsort(d[:, 0])] + if name.startswith("spa") or name == "exact_total": + d = d[:: max(1, len(d) // 80)] # thin the per-column line samples + series.setdefault(name, []).append(d) + return {k: np.vstack(v) for k, v in series.items()} + + +def main(): + PDF.parent.mkdir(exist_ok=True) + if not PDF.exists(): + print(f"downloading {ARXIV_URL} ...") + urllib.request.urlretrieve(ARXIV_URL, PDF) + page = fig6_page(PDF) + png = render(PDF, page, HERE / "data" / "ge_fig6") + series = digitize(png) + + with open(OUT, "w") as f: + f.write("# Ge-Chong-Stone PRA 82, 063824 (arXiv:1008.0628) Fig. 6, digitized\n") + f.write("# exact_* : symbols (numerical solution of Eq. 28, exact SALT)\n") + f.write("# spa_* : solid lines (single-pole approximation, Eqs. 40, 44-45)\n") + f.write("# NB: in the single-mode regime the symbols ride on the lines, so the\n") + f.write("# spa_* series carry small bumps where merged markers were not separable.\n") + f.write("series,D0,intensity\n") + for name in sorted(series): + for D0, inten in series[name]: + f.write(f"{name},{D0:.4f},{inten:.4f}\n") + print(f"wrote {OUT} ({sum(len(v) for v in series.values())} rows)") + + +if __name__ == "__main__": + main() diff --git a/examples/line_fabry_perot/README.md b/examples/line_fabry_perot/README.md new file mode 100644 index 00000000..4fbd9827 --- /dev/null +++ b/examples/line_fabry_perot/README.md @@ -0,0 +1,11 @@ +# Open Fabry–Pérot line + +A 1D path graph whose two end edges are vacuum leads (the radiative loss that +sets the lasing threshold). With the narrow gain used here, two longitudinal +modes near the gain centre lase under both `linear` and `full_salt_newton`; +above threshold the newton curves bend below the linear ones — the full-SALT +gain saturation the near-threshold model omits. + +`bash run.sh` writes `li_curves.png` (geometry | per-mode L–I | total) here. +Figures are not committed; re-run to reproduce. See `doc/source/lasing.rst` +for the solver physics. diff --git a/examples/line_fabry_perot/run.py b/examples/line_fabry_perot/run.py new file mode 100644 index 00000000..44d7c5c9 --- /dev/null +++ b/examples/line_fabry_perot/run.py @@ -0,0 +1,31 @@ +"""Open 1D Fabry--Perot cavity: ``linear`` vs ``full_salt_newton`` L--I curves. + +A path graph whose two end edges are vacuum leads (the radiative loss that sets +the lasing threshold). A short optical length keeps the longitudinal modes well +separated in ``k``; with the narrow gain here two modes near the line centre +lase under both solvers, and the newton curves bend below the linear ones above +threshold (the full-SALT gain saturation the near-threshold model omits). + +Run from this directory (writes ``li_curves.png`` here):: + + OMP_NUM_THREADS=1 python run.py +""" + +import sys +from pathlib import Path + +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import compare_and_plot, quantum_graph + + +def build(n_edges=10, total_length=0.5): + g = nx.path_graph(n_edges + 1) + pos = np.array([[i, 0.0] for i in range(n_edges + 1)], dtype=float) + return quantum_graph(g, pos, total_length) + + +if __name__ == "__main__": + compare_and_plot(build(), "line (Fabry-Perot)", Path(__file__).resolve().parent) diff --git a/examples/line_fabry_perot/run.sh b/examples/line_fabry_perot/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/line_fabry_perot/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/long_line/README.md b/examples/long_line/README.md new file mode 100644 index 00000000..78718348 --- /dev/null +++ b/examples/long_line/README.md @@ -0,0 +1,14 @@ +# Long multimode Fabry–Pérot + +The open 1D cavity of `../line_fabry_perot`, four times longer: the +longitudinal mode spacing dk = pi/(nL) ≈ 0.52 packs ~12 modes into the gain +window and five of them lase — the textbook multimode FP laser, driven by +spatial hole burning (adjacent longitudinal modes have shifted standing-wave +patterns, so each finds gain the others left). `linear` lases 4; +`full_salt_newton` adds a fifth late mode the frozen-profile model suppresses, +while the two solvers' **total** L–I curves coincide — per-mode redistribution +with the total preserved, the same pattern validated against Ge–Chong–Stone on +`../line_PRA`. + +`bash run.sh` writes `li_curves.png` + `mode_profiles.png` here. Figures are +not committed; re-run to reproduce. diff --git a/examples/long_line/run.py b/examples/long_line/run.py new file mode 100644 index 00000000..414a3d26 --- /dev/null +++ b/examples/long_line/run.py @@ -0,0 +1,43 @@ +"""Long multimode Fabry--Perot: many longitudinal modes by cavity length. + +The same open 1D cavity as ``../line_fabry_perot`` but four times longer +(optical length ``nL = 6``): the longitudinal mode spacing +``dk = pi/(nL) ~ 0.52`` packs ~12 modes into the scan window, and the broad +gain (``gamma_perp = 3``) lets several of them reach threshold -- the textbook +multimode Fabry--Perot laser. Adjacent longitudinal modes have shifted +standing-wave patterns, so each burns its spatial holes in different places +and they co-lase despite sharing the cavity (the classic spatial-hole-burning +route to multimode operation in FP diodes). + +Checks encoded here: the mode spacing matches ``pi/(nL)``, and ``linear`` and +``full_salt_newton`` agree on the lasing count near threshold. + +Run from this directory (writes ``li_curves.png`` + ``mode_profiles.png``):: + + OMP_NUM_THREADS=1 python run.py +""" + +import sys +from pathlib import Path + +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import compare_and_plot, quantum_graph + +TOTAL_LENGTH = 2.0 # 4x the short line -> dk = pi/(3*2) ~ 0.52 + + +def build(n_edges=20, total_length=TOTAL_LENGTH): + g = nx.path_graph(n_edges + 1) + pos = np.array([[i, 0.0] for i in range(n_edges + 1)], dtype=float) + return quantum_graph(g, pos, total_length) + + +if __name__ == "__main__": + graph = build() + inner = [graph[u][v]["length"] for u, v in graph.edges if graph[u][v]["inner"]] + dk = np.pi / (3.0 * sum(inner)) + print(f"cavity length {sum(inner):.3f} -> expected longitudinal spacing dk = {dk:.3f}") + compare_and_plot(graph, "long Fabry-Perot line", Path(__file__).resolve().parent) diff --git a/examples/long_line/run.sh b/examples/long_line/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/long_line/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/mini_buffon/README.md b/examples/mini_buffon/README.md new file mode 100644 index 00000000..8ebae9c5 --- /dev/null +++ b/examples/mini_buffon/README.md @@ -0,0 +1,50 @@ +# Mini-buffon network — how to get more modes lasing on a disordered graph + +A shrunk Nat. Commun.-style buffon network (10 random lines, giant component: +39 nodes, 45 edges, 18 radiating lead ends; fixed seed). The spectrum is +genuinely dense — 10 modes in the scan window, ~25 per unit k (Weyl estimate +nL/π ≈ 29) — with disorder-spread losses and a near-degenerate pair at +Δk = 0.006. + +Three measured stages (all encoded in `run.py`): + +| stage | pump | ceiling | linear lases | newton lases | +|---|---|---|---|---| +| low uniform | all 27 inner edges | 0.4 | 2 | 2 | +| high uniform | all 27 inner edges | 1.2 | 5 | 4 | +| shaped | 20 mode-owned edges | 1.2 | 2 | 3 | + +- **Gain clamping suppresses, but not absolutely**: the losers' interacting + thresholds are finite, just 10–40× their bare ones — sweeping the same + uniform pump 3× further buys back 4–5 co-lasing modes. On this network + **pump strength, not pump shaping, is the simplest route to more modes.** +- **Pump shaping selects rather than multiplies**: every shaped pump probed + (greedy low-cross-saturation targets, several target counts/margins) lased + *fewer* modes than uniform at the same ceiling — removing pump area raises + all thresholds faster than the decoupling pays back. Shaping is the + Nat. Commun. lever for choosing *which* mode lases (`netsalt.pump`). +- **Solver cost is not the constraint at this scale**: `full_salt_newton` + takes ~25–40 s per sweep on the ~700-node oversampled work graph (vs ~0.1 s + for `linear`) and tracks `linear`'s sets. +- **This graph hardened the solver**: early runs showed spurious per-mode + kinks (the Δk = 0.004 pair swapping identities at active-set events, the + bootstrap capturing the trivial a = 0 root, an add flipping a stronger + veteran into a weaker newcomer). Fixed in `full_salt_newton` by + linear-slope warm starts, a candidate-spacing cap on the k-window, and + continuity guards keyed on physical invariants (lower-threshold veteran + killed, total-output monotonicity). The coarse and fine grids now agree up + to a genuine **mode crossing** at D0 ≈ 0.6 (~20× threshold) where mode 8 + overtakes mode 6 — there the frozen-field single-pole iteration cannot + resolve the per-mode split (it falls to a spurious lower-total branch, + robustly across relaxation/step size: a method limit). A **total-output + ratchet** enforces the hard physical law (total cannot fall as pump rises) + by holding the collapsing mode — the dominant curve plateaus (flagged) and + the **total L–I stays monotone and physical**. The per-mode magnitudes + across the crossing are at the method's resolution limit; fully resolving + it would need a constant-flux-state SALT solver. + +For many co-lasing modes by *design* (localisation, not pump), see +`../ring_chain`; for why raw spectral density does not help, `../chord_sweep`. + +`bash run.sh` writes the three `*_li_curves.png` / `*_mode_profiles.png` pairs +here (~6 min). Figures are not committed; re-run to reproduce. diff --git a/examples/mini_buffon/run.py b/examples/mini_buffon/run.py new file mode 100644 index 00000000..4b72b965 --- /dev/null +++ b/examples/mini_buffon/run.py @@ -0,0 +1,208 @@ +"""Mini-buffon network: the dense-spectrum stress test for the two solvers. + +A shrunk version of the Nat. Commun. buffon networks (10 random lines, giant +component: 39 nodes, 45 edges, 18 radiating lead ends) built with +:func:`netsalt.utils.make_buffon_graph` from a fixed seed. The spectrum is +genuinely dense -- ~25 modes per unit ``k`` (close to the Weyl estimate +``nL/pi = 28.6``), 10 modes in the scan window with disorder-spread losses +(``alpha = 0.006 - 0.064``) and a near-degenerate pair split by ``dk = 0.006``. + +What the run shows (measured, not assumed) -- three stages: + +* **Low uniform pump (D0 <= 0.4): gain clamping wins.** Ten modes sit under + the broadened gain but only **two** lase -- the extended disorder modes + overlap strongly and the winners clamp the gain. But the suppression is + not absolute: the losers' *interacting* thresholds are finite, just 10-40x + their bare ones. +* **High uniform pump (D0 <= 1.2): pump strength buys the modes back.** The + same uniform pump swept 3x further lases **5 (linear) / 4 (newton)** -- on + this network *more pump*, not pump shaping, is the simplest route to more + co-lasing modes. +* **Pump shaping selects, it does not multiply.** A mode-resolved pump + (greedy low-cross-saturation targets from the competition matrix, pumping + the edges each target dominates) was probed at several target counts and + ownership margins: at the *same* pump ceiling every shaped pump lased + *fewer* modes (2-3) than uniform (4-5), because removing pump area raises + all thresholds faster than the decoupling pays back. Shaping is the tool + for choosing *which* modes lase (the Nat. Commun. optimisation route, + ``netsalt.pump``); the shaped stage is kept here to document that + trade-off honestly. + +Throughout, ``full_salt_newton`` stays comfortable (~25-40 s per sweep on the +~700-node oversampled work graph vs ~0.1 s for ``linear``) -- at this scale +competition physics, not solver cost, is the constraint. For many co-lasing +modes by *design* see ``../ring_chain``. + +**History -- this graph hardened the solver.** Early runs showed spurious +per-mode kinks. Several artifact classes were traced here and fixed in the +solver: the near-degenerate pair (modes 5/6, ``dk = 0.004``) swapping +identities at active-set events; the bootstrap capturing the trivial +``a = 0`` root while its ``k`` drifted onto the twin; an *add* whose +re-solve flipped a stronger veteran's amplitude into a weaker newcomer. The +fixes -- linear-slope warm starts, a candidate-spacing cap on the +``k``-window, and continuity guards keyed on the physical invariants +(lower-threshold-veteran-killed, and total-output monotonicity) -- make the +coarse and fine pump grids agree everywhere up to a genuine **mode +crossing** near ``D0 ~ 0.6`` (~20x threshold), where mode 8 overtakes mode +6. There the frozen-field single-pole iteration cannot resolve the per-mode +split (it converges to a spurious lower-total branch, robustly across +relaxation and step size -- a limit of the method, not a tuning bug). A +**total-output ratchet** enforces the one hard physical law there -- the +total cannot fall as the pump rises -- by holding any collapsing mode at its +last value (the dominant curve plateaus, flagged by a warning) while letting +the others grow. The total L--I is therefore monotone and physical; the +per-mode magnitudes across the crossing are at the method's resolution +limit, and the plateau marks it honestly. Fully resolving the crossing would +need a constant-flux-state SALT solver. + +Run from this directory (writes ``li_curves.png`` + ``mode_profiles.png``):: + + OMP_NUM_THREADS=1 python run.py # ~2 minutes, mostly the mode pipeline +""" + +from __future__ import annotations + +import sys +import time +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import compare_and_plot + +import netsalt +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import make_buffon_graph + +HERE = Path(__file__).resolve().parent +N_LINES = 10 +SEED = 4 # chosen by a small seed scan for a ~40-node giant component +D0_MAX = 0.4 + +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 3.55, # on the low-loss cluster found by the spectrum scan + "gamma_perp": 0.3, + "k_min": 3.3, + "k_max": 3.7, + "alpha_min": -0.05, + "alpha_max": 0.3, + "n_workers": 1, + "n_modes_max": 60, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": D0_MAX, + "D0_steps": 14, + "dielectric_params": {"method": "uniform", "inner_value": 9.0, "outer_value": 1.0, "loss": 0.0}, +} + + +def build(): + g, pos = make_buffon_graph(n_lines=N_LINES, size=(-100.0, 100.0), resolution=100.0, rng=SEED) + giant = max(nx.connected_components(g), key=len) + g = nx.convert_node_labels_to_integers(g.subgraph(giant).copy(), label_attribute="old") + positions = np.array([pos[g.nodes[u]["old"]] for u in g.nodes]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 30.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def selective_pump(graph, tdf, n_targets=4, margin=0.7): + """Mode-resolved pump: cover the edges each low-overlap target dominates. + + Greedy target selection on the (uniform-pump) competition matrix: start + from the lowest-threshold mode, then repeatedly add the candidate whose + worst normalised cross-saturation ``T_ij / sqrt(T_ii T_jj)`` against the + already-selected set is smallest -- spatially distinct modes by + construction. The pump then covers exactly the inner edges whose largest + per-edge intensity (among *all* candidates) belongs to a target: each + target keeps the gain where it is strongest, the non-targets are starved. + """ + from netsalt.modes import compute_mode_competition_matrix, mean_mode_on_edges + + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + tms = tdf["threshold_lasing_modes"].to_numpy() + candidates = [int(i) for i in np.where(thr < np.inf)[0]] + T = compute_mode_competition_matrix(graph, tdf) + + def xsat(i, j): + return abs(T[i, j]) / max(np.sqrt(abs(T[i, i]) * abs(T[j, j])), 1e-12) + + targets = [min(candidates, key=lambda i: thr[i])] + while len(targets) < min(n_targets, len(candidates)): + rest = [c for c in candidates if c not in targets] + targets.append(min(rest, key=lambda c: max(xsat(c, t) for t in targets))) + print(f"selective-pump targets (low mutual cross-saturation): {sorted(targets)}") + + from netsalt.utils import from_complex + + fields = {} + for i in candidates: + e2 = np.abs(mean_mode_on_edges(from_complex(tms[i]), graph, check_quality=False)) + fields[i] = e2 / max(e2.sum(), 1e-12) + inner = np.array([1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()]) + # pump an edge when a *target* is (close to) the strongest candidate on it: + # strict winner-take-all starves the targets of pump area, so allow any + # target within ``margin`` of the top share + top = np.array([max(fields[i][e] for i in candidates) for e in range(len(inner))]) + near = np.array( + [any(fields[t][e] >= margin * top[e] for t in targets) for e in range(len(inner))] + ) + pump = inner * near + print(f"selective pump covers {int(pump.sum())} of {int(inner.sum())} inner edges") + return pump + + +if __name__ == "__main__": + graph = build() + n_leads = sum(1 for n in graph.nodes if len(graph[n]) == 1) + inner = sum(graph[u][v]["length"] for u, v in graph.edges if graph[u][v]["inner"]) + print( + f"mini-buffon: {len(graph)} nodes, {len(graph.edges)} edges, {n_leads} leads; " + f"Weyl density nL/pi = {3 * inner / np.pi:.1f} modes per unit k" + ) + t0 = time.perf_counter() + uniform = compare_and_plot( + graph, + "mini-buffon (uniform pump, low)", + HERE, + d0_max=D0_MAX, + d0_steps=12, + passive_method="contour", + prefix="uniform_low_", + ) + compare_and_plot( + graph, + "mini-buffon (uniform pump, high)", + HERE, + d0_max=3.0 * D0_MAX, + d0_steps=12, + passive_method="contour", + prefix="uniform_high_", + ) + inner_pump = np.array([1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()]) + graph.graph["params"]["pump"] = inner_pump # design against the uniform operator + pump = selective_pump(graph, uniform["tdf"]) + compare_and_plot( + graph, + "mini-buffon (selective pump)", + HERE, + d0_max=3.0 * D0_MAX, + d0_steps=12, + passive_method="contour", + pump=pump, + prefix="selective_", + ) + print(f"total wall time {time.perf_counter() - t0:.0f}s") diff --git a/examples/mini_buffon/run.sh b/examples/mini_buffon/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/mini_buffon/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/ring_chain/README.md b/examples/ring_chain/README.md new file mode 100644 index 00000000..437c895f --- /dev/null +++ b/examples/ring_chain/README.md @@ -0,0 +1,15 @@ +# Detuned ring chain — the mode count as a dial + +Four rings of increasing size bridged in sequence. Detuning localises one mode +per ring near the gain centre (k = 3.61/3.66/3.73/3.83, 87–100 % on the home +ring); the output leads hang off the **bridge midpoints**, so they selectively +damp the hybridised (delocalised) modes while the one-per-ring quartet keeps +only the uniform material-loss floor. Result: **four co-lasing modes, one per +ring** under both `linear` and `full_salt_newton`, with onsets staggered by +loss and gain detuning (ring 2's mode is the leakiest and turns on last). Add +a ring, add a mode. The mode-profile figure shows each mode on its own ring +with near-zero saturated reshaping — weak competition is *why* they co-lase. + +`bash run.sh` writes `ring_chain.png` + `mode_profiles.png` here. Figures are +not committed; re-run to reproduce. See `doc/source/lasing.rst` for the +physics. diff --git a/examples/ring_chain/run.py b/examples/ring_chain/run.py new file mode 100644 index 00000000..5db6f728 --- /dev/null +++ b/examples/ring_chain/run.py @@ -0,0 +1,227 @@ +"""Detuned ring chain: the mode count as a dial (one lasing mode per ring). + +A chain of four rings of increasing size, bridged in sequence. Detuning shifts +each ring's mode comb, so near the gain centre each ring contributes one +**localised** mode (100 % of its intensity on its home ring, found by the +parameter scan below). The design trick is the output coupling: each bridge +carries a midpoint node with a lead. Hybridised (delocalised) modes have weight +on the bridges, localised ring modes avoid them -- so the bridge-mounted leads +*selectively damp the hybrids* (cold-cavity loss 0.02--0.04) while the +one-per-ring quartet keeps only the uniform material-loss floor (~0.004, set by +``dielectric_params['loss']``). With the gain (``k_a = 3.71``, +``gamma_perp = 0.18``) covering the quartet at k = 3.61 / 3.66 / 3.73 / 3.83, +the four localised modes lase one per ring with staggered onsets (ring 2's +mode is the leakiest of the quartet and turns on last), and the +hybrids stay below threshold: **the number of co-lasing modes equals the number +of rings**, by construction. (With leads mounted on the rings instead, the loss +hierarchy inverts -- the hybrids lase first and the story muddies; that failed +design is why the leads sit on the bridges.) + +Both solvers should lase the quartet; ``full_salt_newton`` adds the saturated +mode-profile reshaping (see ``mode_profiles.png``). + +Run from this directory (writes ``ring_chain.png`` + ``mode_profiles.png``):: + + OMP_NUM_THREADS=1 python run.py +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import curves, linear_on_grid, mode_profile_figure + +import netsalt +from netsalt.modes import ( + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + find_passive_modes, + find_threshold_lasing_modes, + mean_mode_on_edges, + pump_trajectories, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import from_complex + +HERE = Path(__file__).resolve().parent +RADII = [0.8, 0.95, 1.1, 1.25] # ~18 % detuning ring to ring +RING_NODES = [7, 8, 9, 10] + +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 3.71, # centred on the one-per-ring quartet (3.61/3.66/3.73/3.83) + "gamma_perp": 0.18, + "k_min": 3.2, + "k_max": 3.9, + "alpha_min": -0.05, + "alpha_max": 0.25, + "n_workers": 1, + "n_modes_max": 40, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 0.08, + "D0_steps": 14, + # the uniform material loss sets the threshold floor of the (otherwise + # essentially lossless) localised ring modes + "dielectric_params": { + "method": "uniform", + "inner_value": 9.0, + "outer_value": 1.0, + "loss": 0.02, + }, +} +D0_MAX = 0.08 +D0_STEPS = 16 + + +def build_ring_chain(): + """Chain of detuned rings, bridges carrying midpoint-mounted leads.""" + g = nx.Graph() + pos = {} + ring_nodes = [] + offset, x0 = 0, 0.0 + for r, n in zip(RADII, RING_NODES, strict=True): + nodes = list(range(offset, offset + n)) + ring_nodes.append(nodes) + g.add_edges_from([(nodes[i], nodes[(i + 1) % n]) for i in range(n)]) + for i, nd in enumerate(nodes): + pos[nd] = [x0 + r * np.cos(2 * np.pi * i / n), r * np.sin(2 * np.pi * i / n)] + x0 += 2.6 * r + offset += n + nid = offset + for a, b in zip(ring_nodes[:-1], ring_nodes[1:], strict=True): + ra = max(a, key=lambda nd: pos[nd][0]) + rb = min(b, key=lambda nd: pos[nd][0]) + mid, lead = nid, nid + 1 + nid += 2 + pos[mid] = [(pos[ra][0] + pos[rb][0]) / 2, (pos[ra][1] + pos[rb][1]) / 2] + pos[lead] = [pos[mid][0], pos[mid][1] - 1.4] + g.add_edge(ra, mid) + g.add_edge(mid, rb) + g.add_edge(mid, lead) + positions = np.array([pos[i] for i in range(len(g))]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 18.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g, ring_nodes + + +def _home_ring(mode, graph, ring_sets): + """Index of the ring holding most of the mode's (length-weighted) intensity.""" + e2 = np.abs(mean_mode_on_edges(mode, graph, check_quality=False)) + weights = [] + for s in ring_sets: + weights.append( + sum( + e2[ei] * graph[u][v]["length"] + for ei, (u, v) in enumerate(graph.edges) + if u in s and v in s + ) + ) + weights = np.array(weights) + frac = weights / max(weights.sum(), 1e-12) + return int(np.argmax(frac)), float(frac.max()) + + +def main(): + graph, ring_nodes = build_ring_chain() + ring_sets = [set(n) for n in ring_nodes] + passive = find_passive_modes(graph, method="contour") + if len(passive) == 0: + raise SystemExit("no passive modes found") + graph.graph["params"]["pump"] = np.array( + [1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()] + ) + trajectories = pump_trajectories(passive, graph, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, graph) + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + tms = tdf["threshold_lasing_modes"].to_numpy() + n_modes = len(tdf) + + homes = {} + print("modes below max pump (id, k, threshold, home ring, localisation):") + for i in range(n_modes): + if thr[i] < D0_MAX: + home, frac = _home_ring(tms[i], graph, ring_sets) + homes[i] = home + print( + f" {i}: k={from_complex(tms[i])[0]:.4f} thr={thr[i]:.4f}" + f" ring {home} ({100 * frac:.0f}%)" + ) + + competition = compute_mode_competition_matrix(graph, tdf) + first = float(thr[thr < np.inf].min()) + grid = np.linspace(first, D0_MAX, D0_STEPS) + linear = linear_on_grid(tdf, competition, grid, n_modes) + n_cols, newton = curves( + compute_modal_intensities_full_salt_newton(graph, tdf.copy(), D0_MAX, D0_steps=D0_STEPS) + ) + peak = max(linear.max(), newton.max(), 1e-9) + ids_linear = [m for m in range(n_modes) if linear[m, -1] > 1e-2 * peak] + ids_newton = [m for m in range(n_modes) if newton[m, -1] > 1e-2 * peak] + rings_lasing = sorted({homes.get(m, -1) for m in ids_newton}) + print(f"linear lases {len(ids_linear)}: {ids_linear}") + print(f"full_salt_newton lases {len(ids_newton)}: {ids_newton} (rings {rings_lasing})") + + cmap = plt.get_cmap("tab10") + fig, axes = plt.subplots(1, 2, figsize=(13, 4.4)) + pos = {n: np.asarray(graph.nodes[n]["position"], dtype=float) for n in graph.nodes} + for u, v in graph.edges(): + x, y = zip(pos[u], pos[v], strict=True) + axes[0].plot(x, y, color="0.55" if graph[u][v]["inner"] else "crimson", lw=2.0) + axes[0].set_aspect("equal") + axes[0].axis("off") + axes[0].set_title(f"{len(RADII)} detuned rings, bridge-mounted leads") + active = sorted(set(ids_linear) | set(ids_newton)) + for m in active: + col = cmap(m % 10) + axes[1].plot(grid, linear[m], "--", color=col, lw=1.3, alpha=0.8) + axes[1].plot( + n_cols, + newton[m], + ".-", + color=col, + lw=1.8, + ms=4, + label=f"mode {m} (ring {homes.get(m, '?')})", + ) + axes[1].set_xlabel("pump $D_0$") + axes[1].set_ylabel("modal intensity") + axes[1].set_title(f"dashed = linear ({len(ids_linear)}), solid = newton ({len(ids_newton)})") + axes[1].legend(fontsize=8) + fig.suptitle("Ring chain: one lasing mode per ring", y=1.02) + fig.tight_layout() + out = HERE / "ring_chain.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") + + mode_profile_figure( + graph, + tdf, + D0_MAX, + ids_linear, + ids_newton, + "ring chain", + HERE, + a0={m: float(newton[m, -1]) for m in range(n_modes)}, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/ring_chain/run.sh b/examples/ring_chain/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/ring_chain/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/ring_leads/README.md b/examples/ring_leads/README.md new file mode 100644 index 00000000..1726c0fd --- /dev/null +++ b/examples/ring_leads/README.md @@ -0,0 +1,10 @@ +# Ring resonator with two leads + +A closed loop made open by two pendant lead edges on opposite sides, so the +loop modes can radiate out. Both `linear` and `full_salt_newton` lase the same +two modes; above threshold the newton curves carry the full-SALT saturation +bend. + +`bash run.sh` writes `li_curves.png` (geometry | per-mode L–I | total) here. +Figures are not committed; re-run to reproduce. See `doc/source/lasing.rst` +for the solver physics. diff --git a/examples/ring_leads/run.py b/examples/ring_leads/run.py new file mode 100644 index 00000000..f547a391 --- /dev/null +++ b/examples/ring_leads/run.py @@ -0,0 +1,34 @@ +"""Ring resonator with two leads: ``linear`` vs ``full_salt_newton`` L--I curves. + +A closed loop made open by two pendant lead edges on opposite sides, so the +loop modes can radiate out. Both solvers lase the same modes; above threshold +the newton curves carry the full-SALT saturation bend. + +Run from this directory (writes ``li_curves.png`` here):: + + OMP_NUM_THREADS=1 python run.py +""" + +import sys +from pathlib import Path + +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import compare_and_plot, quantum_graph + + +def build(n=12, total_length=1.0): + g = nx.cycle_graph(n) + pos = {u: [np.cos(2 * np.pi * u / n), np.sin(2 * np.pi * u / n)] for u in g.nodes} + g.add_edge(0, n) + pos[n] = [2.0, 0.0] + g.add_edge(n // 2, n + 1) + pos[n + 1] = [-2.0, 0.0] + positions = np.array([pos[u] for u in sorted(g.nodes)], dtype=float) + return quantum_graph(g, positions, total_length) + + +if __name__ == "__main__": + compare_and_plot(build(), "ring + leads", Path(__file__).resolve().parent) diff --git a/examples/ring_leads/run.sh b/examples/ring_leads/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/ring_leads/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/tree/README.md b/examples/tree/README.md new file mode 100644 index 00000000..54d75132 --- /dev/null +++ b/examples/tree/README.md @@ -0,0 +1,9 @@ +# Binary-tree splitter + +A depth-3 binary tree: one input lead, several leaf leads. Lossy enough that a +single mode lases under both `linear` and `full_salt_newton` (winner-take-all +gain clamping). + +`bash run.sh` writes `li_curves.png` (geometry | per-mode L–I | total) here. +Figures are not committed; re-run to reproduce. See `doc/source/lasing.rst` +for the solver physics. diff --git a/examples/tree/run.py b/examples/tree/run.py new file mode 100644 index 00000000..bc9e0364 --- /dev/null +++ b/examples/tree/run.py @@ -0,0 +1,29 @@ +"""Binary-tree splitter: ``linear`` vs ``full_salt_newton`` L--I curves. + +A depth-3 binary tree -- one input lead, several leaf leads. Lossy enough that +a single mode lases under both solvers (winner-take-all gain clamping). + +Run from this directory (writes ``li_curves.png`` here):: + + OMP_NUM_THREADS=1 python run.py +""" + +import sys +from pathlib import Path + +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import compare_and_plot, quantum_graph + + +def build(total_length=1.0): + g = nx.balanced_tree(2, 3) + pos = nx.kamada_kawai_layout(g) + positions = np.array([pos[u] for u in sorted(g.nodes)], dtype=float) + return quantum_graph(g, positions, total_length) + + +if __name__ == "__main__": + compare_and_plot(build(), "binary tree", Path(__file__).resolve().parent) diff --git a/examples/tree/run.sh b/examples/tree/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/tree/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/examples/two_ring/README.md b/examples/two_ring/README.md new file mode 100644 index 00000000..cb57afc3 --- /dev/null +++ b/examples/two_ring/README.md @@ -0,0 +1,12 @@ +# Detuned two-ring "photonic molecule" — genuine multimode lasing + +Two rings of different sizes (radii 0.9 / 1.25) joined by a bridge, with a lead +on each. The detuning breaks the left/right symmetry and localises each mode +onto one ring (identical rings would give symmetric/antisymmetric modes spread +over both, with high overlap). Localised modes barely compete, so with a narrow +gain on a cross-ring pair `full_salt_newton` lases several modes spread across +the two rings — genuine multimode lasing, with the modes labelled by ring in +the figure. + +`bash run.sh` writes `two_ring_multimode.png` here. Figures are not committed; +re-run to reproduce. See `doc/source/lasing.rst` for the solver physics. diff --git a/examples/two_ring/run.py b/examples/two_ring/run.py new file mode 100644 index 00000000..24bb9eca --- /dev/null +++ b/examples/two_ring/run.py @@ -0,0 +1,198 @@ +"""Genuine multimode lasing on a detuned two-ring "photonic molecule". + +The simple single-graph examples (``../line_fabry_perot``, ``../ring_leads``, ``../tree``) are all +*single-mode* under faithful SALT: their modes overlap strongly, so the dominant +mode clamps the gain and holds the others below threshold (``full_salt_newton`` +correctly reports one mode). To get +genuine *multimode* lasing you need modes that occupy **different regions** of the +graph, so each burns its own spatial hole and leaves gain for the others. + +This builds two rings of **different sizes** (radii 0.9 and 1.25, i.e. ~1.4x +detuned) joined by a single bridge edge, with a lead on each ring for output +coupling. The detuning is the whole point: + +* If the rings were *identical* the mirror symmetry would force the eigenmodes to + be symmetric/antisymmetric (bonding/antibonding) combinations spread over *both* + rings -- delocalised, high mutual overlap, strong competition. +* Detuning shifts the two rings' mode combs, so at a given frequency only one ring + is near-resonant and the mode **localises** onto that ring. Two localised modes + (one per ring) barely overlap -> they co-lase. + +With a narrow gain centred on a pair of nearby modes from *different* rings, all +four solvers -- including ``full_salt_newton`` -- lase several modes at once. + +Run:: + + OMP_NUM_THREADS=1 python run.py + +Modes are found by Beyn's contour method (robust on this hand-built graph). The +solve is the operator-level Newton; this is the multimode case it was built for. +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _common import mode_profile_figure + +import netsalt +from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + find_passive_modes, + find_threshold_lasing_modes, + mean_mode_on_edges, + pump_trajectories, +) +from netsalt.physics import dispersion_relation_pump +from netsalt.quantum_graph import create_quantum_graph, set_total_length +from netsalt.utils import from_complex + +HERE = Path(__file__).resolve().parent +N_A, N_B = 7, 9 # ring node counts -- different sizes => detuned => localised modes + +# A narrow gain line centred on a cross-ring mode pair near k = 3.57. Mode finding +# uses the contour method, so only k_min/k_max/alpha bounds matter for it. +PARAMS = { + "open_model": "open", + "c": 1.0, + "k_a": 3.567, + "gamma_perp": 0.12, + "k_min": 3.45, + "k_max": 3.66, + "alpha_min": -0.05, + "alpha_max": 0.15, + "n_workers": 1, + "n_modes_max": 10, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 1.0, + "D0_steps": 14, + "dielectric_params": {"method": "uniform", "inner_value": 9.0, "outer_value": 1.0, "loss": 0.0}, +} +D0_MAX = 1.0 + + +def build_two_ring(): + """Two detuned rings joined by a bridge, with a lead on each (open cavity).""" + g = nx.disjoint_union(nx.cycle_graph(N_A), nx.cycle_graph(N_B)) + g.add_edge(0, N_A) # bridge between the rings + g.add_edge(2, N_A + N_B) # lead on ring A + g.add_edge(N_A + 4, N_A + N_B + 1) # lead on ring B + pos = {} + for i in range(N_A): + pos[i] = [-1.6 + 0.9 * np.cos(2 * np.pi * i / N_A), 0.9 * np.sin(2 * np.pi * i / N_A)] + for j in range(N_B): + pos[N_A + j] = [ + 1.7 + 1.25 * np.cos(2 * np.pi * j / N_B), + 1.25 * np.sin(2 * np.pi * j / N_B), + ] + pos[N_A + N_B] = [-1.6, -2.2] + pos[N_A + N_B + 1] = [1.7, -2.6] + positions = np.array([pos[i] for i in range(len(g))]) + create_quantum_graph(g, dict(PARAMS), positions=positions) + set_total_length(g, 9.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + return g + + +def _which_ring(mode, graph): + """'A' or 'B' -- the ring holding most of the mode's intensity.""" + e2 = np.abs(mean_mode_on_edges(mode, graph, check_quality=False)) + ring_a = sum(e2[ei] for ei, (u, v) in enumerate(graph.edges) if u < N_A and v < N_A) + ring_b = sum( + e2[ei] for ei, (u, v) in enumerate(graph.edges) if N_A <= u and N_A <= v < N_A + N_B + ) + return "A" if ring_a > ring_b else "B" + + +def main(): + graph = build_two_ring() + passive = find_passive_modes(graph, method="contour") + if len(passive) == 0: + raise SystemExit("no passive modes found") + graph.graph["params"]["pump"] = np.array( + [1.0 if graph[u][v]["inner"] else 0.0 for u, v in graph.edges()] + ) + trajectories = pump_trajectories(passive, graph, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, graph) + thr = np.asarray(tdf["lasing_thresholds"]).ravel() + threshold_modes = tdf["threshold_lasing_modes"].to_numpy() + rings = {i: _which_ring(threshold_modes[i], graph) for i in range(len(tdf)) if thr[i] < np.inf} + print("lasing modes (id, k, threshold, ring):") + for i, r in rings.items(): + print(f" {i}: k={from_complex(threshold_modes[i])[0]:.3f} thr={thr[i]:.3f} ring {r}") + + competition = compute_mode_competition_matrix(graph, tdf) + solvers = { + "linear": compute_modal_intensities(tdf.copy(), D0_MAX, competition), + "full_salt_newton": compute_modal_intensities_full_salt_newton( + graph, tdf.copy(), D0_MAX, D0_steps=18 + ), + } + cmap = plt.get_cmap("tab10") + fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharex=True) + lasing_at_max = {} + endpoint = {} + for ax, (name, df) in zip(axes.ravel(), solvers.items(), strict=True): + cols = np.array( + sorted(c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities") + ) + data = np.nan_to_num(df[[("modal_intensities", c) for c in cols]].to_numpy(dtype=float)) + peak = max(data.max(), 1e-9) + active = [m for m in range(data.shape[0]) if data[m].max() > 1e-2 * peak] + lasing_at_max[name] = [m for m in active if data[m, -1] > 1e-2 * peak] + endpoint[name] = {m: float(data[m, -1]) for m in active} + for m in active: + ax.plot( + cols, + data[m], + ".-", + color=cmap(m % 10), + label=f"mode {m} (ring {rings.get(m, '?')})", + ) + ax.set_title(f"{name} ({len(active)} lasing)") + ax.set_ylabel("modal intensity") + if active: + ax.legend(fontsize=7) + print( + f"{name}: {len(active)} lasing @max " + + str({int(m): round(float(data[m, -1]), 3) for m in active}) + ) + for ax in axes.ravel(): + ax.set_xlabel("pump $D_0$") + fig.suptitle("Detuned two-ring photonic molecule: genuine multimode lasing", y=1.0) + fig.tight_layout() + out = HERE / "two_ring_multimode.png" + fig.savefig(out, dpi=120, bbox_inches="tight") + plt.close(fig) + print(f"wrote {out}") + + mode_profile_figure( + graph, + tdf, + D0_MAX, + lasing_at_max["linear"], + lasing_at_max["full_salt_newton"], + "two-ring molecule", + HERE, + a0=endpoint["full_salt_newton"], + ) + + +if __name__ == "__main__": + main() diff --git a/examples/two_ring/run.sh b/examples/two_ring/run.sh new file mode 100755 index 00000000..d925ab03 --- /dev/null +++ b/examples/two_ring/run.sh @@ -0,0 +1,5 @@ +#!/bin/bash +export OMP_NUM_THREADS=1 +export NUMEXPR_MAX_THREADS=1 + +python run.py diff --git a/netsalt/modes.py b/netsalt/modes.py index 4d62ffb5..ab72824e 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -32,14 +32,16 @@ find_rough_modes_from_scan, refine_mode, ) -from .physics import gamma, q_value +from .physics import dispersion_relation_pump_saturated, gamma, q_value from .quantum_graph import ( + DENSE_EIG_MAX, construct_incidence_matrix, construct_laplacian, construct_weight_matrix, graph_with_params, graph_with_pump, mode_quality, + oversample_graph, set_wavenumber, ) from .utils import from_complex, get_scan_grid, to_complex @@ -424,30 +426,55 @@ def pump_linear(mode_0, graph, D0_0, D0_1): return from_complex(freq * np.sqrt((1.0 + gamma_overlap * D0_0) / (1.0 + gamma_overlap * D0_1))) -def mode_on_nodes(mode, graph): - """Compute the mode solution on the nodes of the graph.""" +def mode_on_nodes(mode, graph, check_quality=True): + """Compute the mode solution on the nodes of the graph. + + ``check_quality`` (default ``True``) raises if the near-null eigenvalue is + above ``quality_threshold`` -- i.e. ``mode`` is not actually a mode of + ``graph``. The self-consistent / full-SALT solvers evaluate a mode's profile + on a graph pumped *above* that mode's threshold (where the linear operator is + no longer singular at the threshold frequency); they pass + ``check_quality=False`` to take the smallest-eigenvalue field anyway. The + result still reduces continuously to the true mode as the pump returns to + threshold. + """ laplacian = construct_laplacian(to_complex(mode), graph) - min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=1, sigma=0, v0=np.ones(len(graph)), which="LM" - ) + # Dense fast path for small graphs (see DENSE_EIG_MAX): a direct eigensolve is + # several times faster than ARPACK shift-invert at small N and returns the same + # nearest-zero eigenpair (smallest-magnitude eigenvalue and its eigenvector). + # Fall back to ARPACK if the operator is non-finite (a probed ``k`` overflowed), + # since ``np.linalg.eig`` would raise. + dense = laplacian.toarray() if laplacian.shape[0] <= DENSE_EIG_MAX else None + if dense is not None and np.isfinite(dense).all(): + eigenvalues, eigenvectors = np.linalg.eig(dense) + idx = int(np.argmin(np.abs(eigenvalues))) + min_eigenvalue = eigenvalues[idx] + node_solution = eigenvectors[:, idx] + else: + min_eigenvalue_arr, node_solution_arr = sc.sparse.linalg.eigs( + laplacian, k=1, sigma=0, v0=np.ones(len(graph)), which="LM" + ) + min_eigenvalue = min_eigenvalue_arr[0] + node_solution = node_solution_arr[:, 0] + quality_thresh = graph.graph["params"].get("quality_threshold", 1e-4) - if abs(min_eigenvalue[0]) > quality_thresh: + if check_quality and abs(min_eigenvalue) > quality_thresh: raise ValueError( "Not a mode, as quality is too high: " - + str(abs(min_eigenvalue[0])) + + str(abs(min_eigenvalue)) + " > " + str(quality_thresh) + ", mode: " + str(mode) ) - return node_solution[:, 0] + return node_solution -def flux_on_edges(mode, graph): +def flux_on_edges(mode, graph, check_quality=True): """Compute the flux on each edge (in both directions).""" - node_solution = mode_on_nodes(mode, graph) + node_solution = mode_on_nodes(mode, graph, check_quality=check_quality) _, B = construct_incidence_matrix(graph) Winv = construct_weight_matrix(graph, with_k=False) @@ -455,10 +482,18 @@ def flux_on_edges(mode, graph): return Winv.dot(B).dot(node_solution) -def mean_mode_on_edges(mode, graph): +def mean_mode_on_edges(mode, graph, check_quality=True): r"""Compute the average :math:`Real(E^2)` on each edge.""" - edge_flux = flux_on_edges(mode, graph) + edge_flux = flux_on_edges(mode, graph, check_quality=check_quality) + return _mean_intensity_from_flux(edge_flux, graph) + + +def _mean_intensity_from_flux(edge_flux, graph): + """Per-edge average ``|E|^2`` from a precomputed edge flux. + Split out of :func:`mean_mode_on_edges` so callers that already hold the flux + (and ``graph.graph['ks']`` for the same mode) avoid a second eigen-solve. + """ mean_edge_solution = np.zeros(len(graph.edges)) for ei in range(len(graph.edges)): k = 1.0j * graph.graph["ks"][ei] @@ -579,19 +614,19 @@ def compute_gamma_q_values(graph, modes_df, df_entry="passive"): ] -def _precomputations_mode_competition(graph, pump_mask, mode_threshold): +def _precomputations_mode_competition(graph, pump_mask, mode_threshold, check_quality=True): """precompute some quantities for a mode for mode competition matrix""" mode, threshold = mode_threshold graph = graph_with_pump(graph, threshold) - node_solution = mode_on_nodes(mode, graph) + node_solution = mode_on_nodes(mode, graph, check_quality=check_quality) z_matrix = compute_z_matrix(graph) BT, Bout = construct_incidence_matrix(graph) Winv = construct_weight_matrix(graph, with_k=False) pump_norm = _graph_norm(BT, Bout, Winv, z_matrix, node_solution, pump_mask) - edge_flux = flux_on_edges(mode, graph) / np.sqrt(pump_norm) + edge_flux = flux_on_edges(mode, graph, check_quality=check_quality) / np.sqrt(pump_norm) k_mu = graph.graph["ks"] gam = gamma(to_complex(mode), graph.graph["params"]) @@ -679,30 +714,38 @@ def _compute_mode_competition_element(lengths, params, data, with_gamma=True): return matrix_element -def compute_mode_competition_matrix(graph, modes_df, with_gamma=True): - """Compute the mode competition matrix, or T matrix.""" - threshold_modes = modes_df["threshold_lasing_modes"].to_numpy() - lasing_thresholds_all = modes_df["lasing_thresholds"].to_numpy() - - threshold_modes = threshold_modes[lasing_thresholds_all < np.inf] - lasing_thresholds = lasing_thresholds_all[lasing_thresholds_all < np.inf] - +def _mode_competition_matrix_block( + graph, threshold_modes, pumps, with_gamma=True, check_quality=True +): + """Build the dense competition matrix over the lasing modes only. + + ``pumps`` is the per-mode pump strength at which each mode's profile is + evaluated. The production matrix uses each mode's own lasing threshold + (approximation #2); the self-consistent path passes a single operating pump + for every mode. ``_precomputations_mode_competition`` already treats its + ``threshold`` argument as the pump (``graph_with_pump``), so no change to the + precompute kernel is needed — only the pump values fed in. ``check_quality`` + is forwarded to ``mode_on_nodes`` (set ``False`` when evaluating profiles at + a pump above the modes' thresholds). + """ precomp = partial( _precomputations_mode_competition, graph, _get_mask_matrices(graph.graph["params"])[1], + check_quality=check_quality, ) - chunksize = max(1, int(0.1 * len(lasing_thresholds) / graph.graph["params"]["n_workers"])) - with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool: + n_workers = graph.graph["params"]["n_workers"] + chunksize = max(1, int(0.1 * len(pumps) / n_workers)) + with multiprocessing.Pool(n_workers) as pool: precomp_results = list( tqdm( pool.imap( precomp, - zip(threshold_modes, lasing_thresholds, strict=True), + zip(threshold_modes, pumps, strict=True), chunksize=chunksize, ), - total=len(lasing_thresholds), + total=len(pumps), ) ) @@ -717,8 +760,8 @@ def compute_mode_competition_matrix(graph, modes_df, with_gamma=True): ] ) - chunksize = max(1, int(0.1 * len(input_data) / graph.graph["params"]["n_workers"])) - with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool: + chunksize = max(1, int(0.1 * len(input_data) / n_workers)) + with multiprocessing.Pool(n_workers) as pool: output_data = list( tqdm( pool.imap( @@ -744,18 +787,33 @@ def compute_mode_competition_matrix(graph, modes_df, with_gamma=True): mode_competition_matrix[mu, nu] = output_data[index] index += 1 - pool.close() + return np.real(mode_competition_matrix) - mode_competition_matrix_full = np.zeros( - [ - len(modes_df["threshold_lasing_modes"]), - len(modes_df["threshold_lasing_modes"]), - ] + +def _scatter_competition_block(block, lasing_mask, n_total): + """Place a lasing-only competition block back into a full n_total matrix.""" + full = np.zeros([n_total, n_total]) + full[np.ix_(lasing_mask, lasing_mask)] = block + return full + + +def compute_mode_competition_matrix(graph, modes_df, with_gamma=True): + """Compute the mode competition matrix, or T matrix. + + Each mode's profile is evaluated at its own lasing threshold pump (the + linearised, near-threshold model). + """ + threshold_modes_all = modes_df["threshold_lasing_modes"].to_numpy() + lasing_thresholds_all = modes_df["lasing_thresholds"].to_numpy() + lasing_mask = lasing_thresholds_all < np.inf + + threshold_modes = threshold_modes_all[lasing_mask] + lasing_thresholds = lasing_thresholds_all[lasing_mask] + + block = _mode_competition_matrix_block( + graph, threshold_modes, lasing_thresholds, with_gamma=with_gamma ) - mode_competition_matrix_full[ - np.ix_(lasing_thresholds_all < np.inf, lasing_thresholds_all < np.inf) - ] = np.real(mode_competition_matrix) - return mode_competition_matrix_full + return _scatter_competition_block(block, lasing_mask, len(threshold_modes_all)) def _find_next_lasing_mode( @@ -793,8 +851,27 @@ def _find_next_lasing_mode( return next_lasing_mode_id, next_lasing_threshold +def _intensity_slopes_shifts(mode_competition_matrix, lasing_thresholds, lasing_mode_ids): + """Linear modal-intensity solve for the active set. + + Returns ``(slopes, shifts)`` such that the modal intensities at pump ``D0`` + are ``slopes * D0 - shifts``. + """ + mode_competition_matrix_inv = np.linalg.pinv( + mode_competition_matrix[np.ix_(lasing_mode_ids, lasing_mode_ids)] + ) + slopes = mode_competition_matrix_inv.dot(1.0 / lasing_thresholds[lasing_mode_ids]) + shifts = mode_competition_matrix_inv.sum(1) + return slopes, shifts + + def compute_modal_intensities(modes_df, max_pump_intensity, mode_competition_matrix): - """Compute the modal intensities of the modes up to D0, with D0_steps.""" + """Compute the modal intensities of the modes up to D0, with D0_steps. + + Event-driven sweep over the pump strength with the fixed (near-threshold) + competition matrix: intensities grow piecewise-linearly between mode + activation / vanishing events. + """ lasing_thresholds = np.asarray(modes_df["lasing_thresholds"]).ravel() next_lasing_mode_id = int(np.argmin(lasing_thresholds)) @@ -810,25 +887,35 @@ def compute_modal_intensities(modes_df, max_pump_intensity, mode_competition_mat pump_intensity = next_lasing_threshold L.debug("Max pump intensity %s", max_pump_intensity) + # safety cap so the event loop always terminates (it needs <~2*n_modes events) + max_events = 100 * (len(modes_df) + 1) + event = 0 while pump_intensity <= max_pump_intensity: + event += 1 + if event > max_events: + warnings.warn( + "modal-intensity sweep hit its event cap; returning the partial L--I curve.", + stacklevel=2, + ) + break L.debug("Current pump intensity %s", pump_intensity) # 1) compute the current mode intensities - mode_competition_matrix_inv = np.linalg.pinv( - mode_competition_matrix[np.ix_(lasing_mode_ids, lasing_mode_ids)] + slopes, shifts = _intensity_slopes_shifts( + mode_competition_matrix, lasing_thresholds, lasing_mode_ids ) - slopes = mode_competition_matrix_inv.dot(1.0 / lasing_thresholds[lasing_mode_ids]) - shifts = mode_competition_matrix_inv.sum(1) # if we hit the max intensity, we add last points and stop if pump_intensity >= max_pump_intensity: L.debug("Max pump intensity reached.") - modal_intensities.loc[lasing_mode_ids, max_pump_intensity] = ( - slopes * max_pump_intensity - shifts + modal_intensities.loc[lasing_mode_ids, max_pump_intensity] = np.clip( + slopes * max_pump_intensity - shifts, 0.0, None ) break - modal_intensities.loc[lasing_mode_ids, pump_intensity] = slopes * pump_intensity - shifts + modal_intensities.loc[lasing_mode_ids, pump_intensity] = np.clip( + slopes * pump_intensity - shifts, 0.0, None + ) # 2) search for next lasing mode next_lasing_mode_id, next_lasing_threshold = _find_next_lasing_mode( @@ -894,6 +981,818 @@ def compute_modal_intensities(modes_df, max_pump_intensity, mode_competition_mat return modes_df +def _finalise_modal_intensities(modes_df, modal_intensities, interacting_lasing_thresholds): + """Attach an L--I sweep to ``modes_df`` (shared by the iterative solvers). + + Mirrors the tail of :func:`compute_modal_intensities`: stores the + interacting thresholds and one ``("modal_intensities", D0)`` column per pump, + rounded to 8 decimals. Columns are written in increasing pump order so the + L--I curve is monotone in ``D0`` regardless of the order they were computed. + """ + modes_df["interacting_lasing_thresholds"] = interacting_lasing_thresholds + + if "modal_intensities" in modes_df: + del modes_df["modal_intensities"] + + pumps = sorted(modal_intensities.columns) + for pump_intensity in pumps: + modes_df["modal_intensities", np.around(pump_intensity, 8)] = modal_intensities[ + pump_intensity + ] + + n_lasing = 0 + if pumps: + last = np.nan_to_num(modal_intensities[pumps[-1]].to_numpy()) + n_lasing = int(np.sum(last > 0)) + L.info("%s lasing modes out of %s", n_lasing, len(modal_intensities.index)) + return modes_df + + +def _single_mode_field_intensity(graph, mode, pump_mask): + """Normalised per-edge field intensity of ``mode`` on ``graph``. + + Returns ``mean_mode_on_edges`` (per-edge ``|E|^2``) divided by the + pump-region norm ``∫_pump |ψ|^2`` (``_graph_norm``) -- the same normalisation + the linear competition matrix applies to its edge fluxes + (``_precomputations_mode_competition``). The amplitude ``a`` it scales is thus + a modal-intensity in the ``∫_pump |Ê|^2 = 1`` convention (its own unit; see + :func:`compute_modal_intensities_full_salt_newton`). + """ + # Nudge alpha off exactly 0: at the lasing point the operator is singular and + # ARPACK shift-invert (sigma=0) cannot factorise it; a tiny imaginary offset + # lifts the singularity and the field is continuous there. + mode = [float(mode[0]), float(mode[1]) if abs(mode[1]) > 1e-7 else 1e-7] + node_solution = mode_on_nodes(mode, graph, check_quality=False) + z_matrix = compute_z_matrix(graph) + BT, Bout = construct_incidence_matrix(graph) + Winv = construct_weight_matrix(graph, with_k=False) + pump_norm = _graph_norm(BT, Bout, Winv, z_matrix, node_solution, pump_mask) + # reuse the single eigen-solve above for the per-edge intensity (avoids a + # second mode_on_nodes inside mean_mode_on_edges) + edge_flux = Winv.dot(Bout).dot(node_solution) + intensity = _mean_intensity_from_flux(edge_flux, graph) + return intensity / abs(pump_norm) + + +def _saturated_graph_at(graph, mode, a, D0, pump, field_intensity): + """Throwaway copy whose ``D0_eff`` carries the spatial-hole-burning denominator. + + ``D0_eff[e] = D0·pump[e] / (1 + Γ(k)·a·|Ê(x)|^2[e])`` with the Lorentzian + gain clamp ``Γ(k) = -Im γ(k)``. The copy uses + :func:`dispersion_relation_pump_saturated`, so ``construct_laplacian`` builds + the saturated operator ``L_sat(k)`` from it without touching shared state. + """ + gain_clamp = -np.imag(gamma(to_complex(mode), graph.graph["params"])) + denom = 1.0 + gain_clamp * a * np.asarray(field_intensity) + g = graph_with_params(graph, D0_eff=D0 * np.asarray(pump) / denom) + g.graph["dispersion_relation"] = dispersion_relation_pump_saturated + return g + + +def _refine_local(mode, graph, tol, max_steps, seed, k_window=1.0): + """Local complex-``k`` refine driving ``(Re λ₁, Im λ₁) → 0`` via MINPACK ``hybr``. + + The engine of *continuous mode-following*: started from the mode's previous + position it converges to the nearby root, so each mode tracks itself across + pump and ARPACK cannot swap one mode's eigenvalue for another's. A fixed + ``default_rng(seed)`` makes the ARPACK start vector (hence the residual) + deterministic; an exactly-singular factorisation (the structural ``k = 0`` DC + point, or the lasing point itself) is read a hair off in ``k``; and a refined + mode that jumped more than ``k_window`` is rejected (it left its branch). + """ + mode = np.asarray(mode, dtype=float) + + def residual(x): + try: + lam = mode_quality( + x, graph, quality_method="complex_eigenvalue", rng=np.random.default_rng(seed) + ) + except RuntimeError: + lam = mode_quality( + [x[0] + 1e-7, x[1]], + graph, + quality_method="complex_eigenvalue", + rng=np.random.default_rng(seed), + ) + return [lam.real, lam.imag] + + result = sc.optimize.root( + residual, + mode, + method="hybr", + tol=0, + options={"maxfev": int(max_steps), "xtol": max(tol, 1e-9)}, + ) + refined = np.asarray(result.x) + if abs(refined[0] - mode[0]) > k_window: + return mode + return refined + + +def _d0_eff_array(graph, modes, a_arr, D0, pump, fields): + """Per-edge saturated effective pump ``D0_eff`` for the multi-mode field.""" + pump = np.asarray(pump) + denom = np.ones(len(pump)) + for m, a_nu, f_nu in zip(modes, a_arr, fields, strict=True): + gain_clamp = -np.imag(gamma(to_complex(m), graph.graph["params"])) + denom = denom + gain_clamp * a_nu * np.asarray(f_nu) + return D0 * pump / denom + + +def _saturated_graph_from_d0_eff(graph, d0_eff): + """Throwaway copy carrying a precomputed ``D0_eff`` and the saturated law.""" + g = graph_with_params(graph, D0_eff=d0_eff) + g.graph["dispersion_relation"] = dispersion_relation_pump_saturated + return g + + +def _saturated_graph_multi(graph, modes, a_arr, D0, pump, fields): + """Throwaway copy carrying the *multi-mode* spatial-hole-burning denominator. + + ``D0_eff[e] = D0·pump[e] / (1 + Σ_ν Γ(k_ν)·a_ν·|Ê_ν(x)|^2[e])`` -- every + lasing mode burns the shared gain. Reduces to :func:`_saturated_graph_at` for + one mode. + """ + return _saturated_graph_from_d0_eff(graph, _d0_eff_array(graph, modes, a_arr, D0, pump, fields)) + + +def _lam_real_k(graph, k, seed): + """Complex ``λ₁`` of ``graph`` at the real wavenumber ``k``. + + Evaluated only on the real axis, where the saturated lasing operator is + (near-)singular and ARPACK-friendly; an exactly-singular factorisation is read + a hair off in ``k``. + """ + try: + return mode_quality( + [k, 0.0], graph, quality_method="complex_eigenvalue", rng=np.random.default_rng(seed) + ) + except RuntimeError: + return mode_quality( + [k + 1e-7, 0.0], + graph, + quality_method="complex_eigenvalue", + rng=np.random.default_rng(seed), + ) + + +def _salt_block_residual(graph, x, n, D0, pump, fields, seed): + """Stacked ``[Re λ₁, Im λ₁]`` at real ``k_μ`` for the active set, fields frozen. + + Holding the hole-burning fields fixed makes each residual a *single* eigensolve + per mode (no inner fixed point) -- a clean, noise-free function of ``(k, a)`` + whose finite-difference Jacobian the trust region can rely on. + """ + ks = x[:n] + a = np.clip(x[n:], 0.0, None) + g = _saturated_graph_multi(graph, [[k, 0.0] for k in ks], a, D0, pump, fields) + out = np.empty(2 * n) + for i in range(n): + lam = _lam_real_k(g, ks[i], seed) + out[2 * i], out[2 * i + 1] = lam.real, lam.imag + return out + + +def _solve_active_set( + graph, + modes0, + fields0, + a0, + D0, + pump, + pump_mask, + max_steps, + seed, + outer=6, + damping=0.7, + k_window_cap=None, +): + """Frozen-field trust-region ``(k, a)`` solve for a *fixed* active set. + + Block iteration: (i) freeze the saturated background fields, (ii) solve every + mode's ``(k_μ real, a_μ ≥ 0)`` with one bounded trust-region least-squares on + the clean :func:`_salt_block_residual` (``k`` confined to a window below the + inter-mode spacing so modes cannot hop, ``a`` to ``[0, a_max]``), (iii) refresh + the fields, repeat. This replaces the old decoupled amplitude least-squares + whose residual re-ran an inner fixed point -- a noisy Jacobian that made the + multimode amplitudes chatter. Returns ``(modes, fields, a, converged)``. + """ + n = len(modes0) + k0 = np.array([float(m[0]) for m in modes0]) + a = np.clip(np.asarray(a0, dtype=float), 1e-3, None) + fields = [np.asarray(f, dtype=float) for f in fields0] + # Confine k to a *tight* window: frequency pulling above threshold is small, + # and a loose window lets the trust region zero a mode's residual by drifting + # its k to a spurious nearby root with a=0 (collapsing multimode to one mode) + # instead of raising its amplitude to lase. The window tracks the spacing + # (0.2 * min gap) with only a numerical floor, so it stays below the spacing + # even for a dense/near-degenerate spectrum (a fixed floor that exceeded the + # spacing made near-degenerate modes collide -> collapse or divergence). + if n > 1: + gaps = np.abs(k0[:, None] - k0[None, :]) + gaps[np.diag_indices(n)] = np.inf + window = float(np.clip(0.2 * gaps.min(), 1e-6, 0.1)) + else: + window = 0.1 + if k_window_cap is not None: + # on a dense spectrum a single active mode's fixed 0.1 window spans many + # roots; cap every window by the spacing of the *full* candidate set so + # a solve cannot drift a mode onto a neighbour's root (observed on the + # mini-buffon near-degenerate pair: the bootstrap solve captured the + # trivial a = 0 branch while k slid onto the twin) + window = min(window, float(k_window_cap)) + a_max = max(1.0e3 * max(float(np.max(a)), 1.0e-3), 1.0e3) + lo = np.concatenate([k0 - window, np.zeros(n)]) + hi = np.concatenate([k0 + window, np.full(n, a_max)]) + ks = k0.copy() + converged = False + for _ in range(outer): + x0 = np.clip(np.concatenate([ks, a]), lo, hi) + try: + result = sc.optimize.least_squares( + lambda x, _f=fields: _salt_block_residual(graph, x, n, D0, pump, _f, seed), + x0, + bounds=(lo, hi), + xtol=1e-6, + ftol=1e-6, + gtol=1e-6, + max_nfev=int(max_steps), + ) + ks, a = result.x[:n], np.clip(result.x[n:], 0.0, None) + except (RuntimeError, ValueError, sc.sparse.linalg.ArpackError): + break + g = _saturated_graph_multi(graph, [[k, 0.0] for k in ks], a, D0, pump, fields) + new = [_single_mode_field_intensity(g, [ks[i], 0.0], pump_mask) for i in range(n)] + change = sum(np.linalg.norm(new[i] - fields[i]) for i in range(n)) + fields = [(1.0 - damping) * fields[i] + damping * new[i] for i in range(n)] + if change <= 1e-6 * (1.0 + sum(np.linalg.norm(f) for f in fields)): + converged = True + break + modes = [np.array([float(ks[i]), 0.0]) for i in range(n)] + return modes, fields, a, converged + + +def _newton_onset_unit_scale(graph, mode0, field0, threshold, t_self): + """Per-mode factor converting the Newton amplitude to the linear-intensity unit. + + The Newton amplitude solves the *operator* clamp; the linear/SPA intensity uses + the competition diagonal ``T_μμ``. They share the onset slope once the Newton + amplitude is rescaled by ``s_linear / s_newton``, with ``s_linear = + 1/(T_μμ·D0_thr)`` and ``s_newton = da/dD0`` the Newton amplitude's slope *at* + threshold. ``s_newton`` is **analytic** -- first-order perturbation of the + lasing condition. The saturated operator perturbs the effective gain + ``g = γ·(D0·P - D0·Γ·a·H)`` with the coherent overlap factors + + .. math:: P = \\frac{\\int_{pump} E^2}{\\int_{inner} ε E^2}, \\qquad + H = \\frac{\\int_{pump} f\\,E^2}{\\int_{inner} ε E^2}, + + where ``f`` is the per-edge hole-burning intensity ``field0`` (per-edge + *constant* on the oversampled work graph -- exactly the profile the operator + itself clamps with, so the first order is exact, not a within-edge + approximation). The mode frequency responds as ``δω ∝ -δg/(1+g)`` + (:func:`pump_linear`); holding ``Im ω = 0`` (the lasing condition) gives + + .. math:: s_{newton} = \\frac{{\\rm Im}[γP/Q]}{D_0^{thr}\\,Γ\\,{\\rm Im}[γH/Q]}, + \\qquad Q = 1 + γ D_0^{thr} P . + + With the shared ``∫_pump |E|^2 = 1`` normalization this comes out ≈ 1 -- + the Newton amplitude *is* the linear modal intensity to first order -- so the + rescale is a small consistency correction, and the near-threshold agreement + between the two solvers is a genuine prediction rather than a calibration + (the earlier *measured* probe at ``1.2·D0_thr`` carried a secant bias of up + to a few % on modes whose curve already bends there). + """ + s_linear = 1.0 / (t_self * threshold) if (t_self > 0.0 and threshold > 0.0) else 0.0 + if s_linear <= 0.0: + return 1.0 + g = graph_with_pump(graph, threshold) + params = g.graph["params"] + node_solution = mode_on_nodes(mode0, g, check_quality=False) + z_matrix = compute_z_matrix(g) + BT, Bout = construct_incidence_matrix(g) + Winv = construct_weight_matrix(g, with_k=False) + inner = np.asarray(params["inner"], dtype=float) + eps_mask = _get_dielectric_constant_matrix(params).dot(sc.sparse.diags(_convert_edges(inner))) + pump_profile = np.asarray(params["pump"], dtype=float) * inner + hole_profile = pump_profile * np.asarray(field0, dtype=float) + inner_norm = _graph_norm(BT, Bout, Winv, z_matrix, node_solution, eps_mask) + p_overlap = ( + _graph_norm( + BT, Bout, Winv, z_matrix, node_solution, sc.sparse.diags(_convert_edges(pump_profile)) + ) + / inner_norm + ) + h_overlap = ( + _graph_norm( + BT, Bout, Winv, z_matrix, node_solution, sc.sparse.diags(_convert_edges(hole_profile)) + ) + / inner_norm + ) + gam = gamma(to_complex(mode0), params) + gain_clamp = -np.imag(gam) + q_factor = 1.0 + gam * threshold * p_overlap + denom = threshold * gain_clamp * np.imag(gam * h_overlap / q_factor) + if not (np.isfinite(denom) and abs(denom) > 0.0): + return 1.0 + s_newton = float(np.imag(gam * p_overlap / q_factor) / denom) + if not (np.isfinite(s_newton) and s_newton > 0.0): + return 1.0 + return float(s_linear / s_newton) + + +NEWTON_DENSE_EIG_MAX = ( + 50 # full_salt_newton runs its (banded, oversampled) eigensolves through ARPACK +) + + +def _auto_oversample_size(graph, modes_df, resolution=12, node_cap=3000): + """Sub-edge length that resolves the lasing standing wave (for hole burning). + + The operator-level hole burning samples ``|E_ν(x)|^2`` per edge; with the bare + edges (one sample per edge) the per-edge **mean** intensity over-estimates the + spatial overlap between modes -- it washes out the within-edge nodes/antinodes + where the coherent competition is actually weak -- so the saturation + **over-clamps** and suppresses co-lasing modes that the competition matrix (and + the Ge-Chong-Stone single-pole SALT, PRA 82, 063824) correctly lases. Sampling + a few points per wavelength fixes it. The local wavelength is + ``λ = 2π / (n·Re k)`` with ``n = sqrt(ε)``; target ``λ_min / resolution``, + capped so the oversampled graph stays bounded. ``resolution`` defaults to 12 + (~λ/12): a convergence study on ``line_PRA`` shows the modal *intensities* + converge to ~1% there, while the cheaper λ/6 (used before the ARPACK + eigensolve scaling) got the lasing count right but was ~15% under-resolved. + ARPACK makes λ/12 essentially free. + """ + cand = np.where(np.asarray(modes_df["lasing_thresholds"]).ravel() < np.inf)[0] + tms = modes_df["threshold_lasing_modes"].to_numpy() + k_max = max((abs(from_complex(tms[i])[0]) for i in cand), default=0.0) + if k_max <= 0.0: + return None + eps = [abs(graph[u][v].get("dielectric_constant", 1.0) or 1.0) for u, v in graph.edges] + n_max = float(np.sqrt(max(eps) if eps else 1.0)) + target = 2.0 * np.pi / (n_max * k_max) / resolution + lengths = np.array([graph[u][v]["length"] for u, v in graph.edges], dtype=float) + est_nodes = float(np.sum(np.maximum(lengths / max(target, 1e-12), 1.0))) + if est_nodes > node_cap: + # Keep the oversampled graph (and its eigensolves) bounded. On a large + # graph this *reduces* the effective resolution below ``resolution`` -- + # raise ``node_cap`` (the full_salt_newton ``oversample_node_cap`` knob) + # to recover within-edge accuracy at higher cost. + target *= est_nodes / node_cap + return float(target) + + +def compute_modal_intensities_full_salt_newton(*args, **kwargs): + """Operator-level full-SALT L--I curves (public entry). + + Thin wrapper that temporarily lowers ``DENSE_EIG_MAX`` so the saturated + eigensolves -- on the *banded*, oversampled graph, targeting isolated lasing + modes -- run through ARPACK shift-invert (~flat in N, far cheaper than dense + O(N^3) at the medium/large sizes oversampling produces). The global default is + left high so dense-spectrum mode *finding* on 2D graphs keeps the robust dense + path. See :data:`~netsalt.quantum_graph.DENSE_EIG_MAX`. + """ + from netsalt import quantum_graph as _qg + + saved = _qg.DENSE_EIG_MAX + _qg.DENSE_EIG_MAX = min(saved, NEWTON_DENSE_EIG_MAX) + try: + return _full_salt_newton_impl(*args, **kwargs) + finally: + _qg.DENSE_EIG_MAX = saved + + +def _full_salt_newton_impl( + graph, + modes_df, + max_pump_intensity, + D0_steps=30, + max_iter=30, + tol=1e-8, + oversample_size=None, + inner_max_iter=25, + inner_damping=0.8, + seed=42, + quality_method="eigenvalue", + oversample_resolution=12, + oversample_node_cap=3000, +): + r"""Operator-level full-SALT L--I curves with a self-consistent active set. + + Solves the real nonlinear SALT eigenproblem: at each pump it finds, for every + *lasing* mode, the real frequency ``k_μ`` and amplitude ``a_μ ≥ 0`` such that + the shared saturated operator ``L_sat`` + (:func:`~netsalt.physics.dispersion_relation_pump_saturated`) is singular at + each real ``k_μ`` simultaneously. + + Two ingredients make the multimode solve robust (see ``doc/source/lasing.rst``): + + * **Frozen-field trust-region solve** (:func:`_solve_active_set`). For a fixed + active set the background fields are frozen while a bounded trust-region + least-squares solves all ``(k_μ, a_μ)``; the fields are then refreshed and the + step repeated. The frozen field makes the residual a single clean eigensolve + per mode, so the Jacobian is noise-free -- unlike the old decoupled solve + whose residual re-ran an inner fixed point and chattered. + * **Self-consistent active set.** The pump is stepped up; at each step the + confirmed lasing set is solved, modes whose amplitude vanishes are dropped, + and non-lasing candidates are probed on the *current saturated background*: + one is added when it has net gain there (``α < 0`` -- the crossing of its + *interacting* threshold), **one per sweep, most above-threshold first** + (simultaneous adds hand the coupled solve a multistable warm start that can + converge to the wrong basin and zero the true winner). Each add is vetted + by a **continuity guard** against identity theft: for near-degenerate + pairs (``dk`` below the gain linewidth) the amplitude split is + ill-conditioned and a *weaker* newcomer can silently steal a *stronger* + veteran's amplitude, which reads as a spurious kink in the L--I curve. + The guard fires when a higher-threshold newcomer kills (``a <= 1e-3``) a + lower-threshold veteran outright -- a genuine takeover leaves the veteran + alive at reduced amplitude, or is led by the lower-threshold mode -- and + then reverts the veterans, rejects the newcomer at this pump step + (retried at the next), and re-probes. Pump-step transitions get the same + continuity treatment with bisected sub-stepping + (``_advance_active_set``). An empty set is + bootstrapped from the lowest noninteracting threshold directly, which is + exact on the unsaturated background (the gain probe alone would reject a + mode sitting *at* its threshold, where ``α = 0``). The active set is thus + found self-consistently from the saturated operator, not borrowed from the + linear model. (This is only faithful when the hole burning is resolved -- + see below; with the bare-edge mean it over-clamps and drops modes that + should co-lase.) + + Amplitudes are reported in the **linear modal-intensity unit** via + :func:`_newton_onset_unit_scale`, the *analytic* (first-order perturbation) + onset slope of the saturated operator. The scale comes out ≈ 1 -- the Newton + amplitude is the linear modal intensity to first order -- so the + near-threshold agreement with ``linear`` is a genuine prediction, not a + calibration. + + **Relation to the SPA competition-matrix solver.** This is the *operator-level* + (exact-spatial) SALT: it solves the real nonlinear eigenproblem rather than the + single-pole-approximation matrix equation + ``D0/D0_thr - 1 = Σ_ν Γ_ν χ_μν I_ν`` (Ge-Chong-Stone, PRA 82, 063824) that + ``linear`` implements. It contributes the + self-consistent gain-clamping active set and the lasing frequencies ``k_μ`` from + the saturated operator, and above threshold it gives the + genuine full-SALT correction beyond the SPA: the dominant mode picks up a + **negative kink** (suppressed *below* the SPA when a second mode turns on and + steals gain) while the second mode sits **above** the SPA, the two nearly + cancelling in the total. Validated against the exact-SALT data of Ge-Chong-Stone + Fig. 6 on ``line_PRA``: the per-mode intensities track the digitized exact curves + to a few percent (dominant 0.21 vs 0.205, second 0.10 vs 0.108 at + ``D0 = 1.27``), and the single-mode regime reduces to the SPA (slope ratio ≈ 1; + see ``examples/line_PRA/compare_to_pra_fig6.py``). + + **Within-edge hole burning must be resolved.** The saturation samples + ``|E_ν(x)|^2`` per edge; with one sample per edge the per-edge *mean* + over-estimates the spatial overlap (it washes out the standing-wave + nodes/antinodes) and **over-clamps**, spuriously suppressing co-lasing modes -- + it lased one mode on ``line_PRA`` where Ge-Chong-Stone (PRA 82, 063824, Eq. 28) + and the competition matrix lase two. ``oversample_size=None`` therefore + auto-picks a wavelength-resolving sub-edge size (:func:`_auto_oversample_size`); + with it newton reproduces the two-mode result and reduces to linear near + threshold. Pass ``oversample_size=0`` for the old (over-clamping) bare-edge + behaviour, or a float to set it explicitly. + + **Scales to buffon networks.** The cost is the oversampled eigensolve, which + ARPACK keeps ~flat in the node count on the banded laplacian, and the + oversampling is bounded by ``oversample_node_cap`` -- so the operator stays at + a few thousand nodes regardless of the cavity's physical length. The real + buffon (96-edge network, ``inner_total_length = 2500``) runs in ~12 s at the + default cap and ~6 min uncapped at ``λ/4`` (≈30k nodes), lasing its + co-lasing modes. Raising ``oversample_node_cap`` (and/or + ``oversample_resolution``) trades speed for within-edge accuracy; on a large + graph the default cap reduces the effective resolution below + ``oversample_resolution``, so pass a higher cap when the modal magnitudes + matter. Two limits to keep in mind on a genuinely *dense* spectrum: (i) + **cost** grows with the number of *co-lasing* modes (the coupled ``(k, a)`` + Jacobian), so it suits narrow gain / pump-targeted few-mode operation rather + than hundreds of simultaneous modes; (ii) **near-degeneracy** -- when the + spacing is far below the gain linewidth the per-mode amplitudes become + ill-conditioned (the total-output ratchet keeps the L--I monotone there). The + ``linear`` competition-matrix solver remains the cheap first pass for the + lasing count at any scale. + + It never raises -- a step that fails to fully converge keeps its iterate and + warns. (``max_iter``, ``tol``, ``inner_max_iter``, ``inner_damping`` are + accepted for interface parity.) + """ + del max_iter, tol, inner_max_iter, inner_damping, quality_method # interface parity + + lasing_thresholds = np.asarray(modes_df["lasing_thresholds"]).ravel() + n_modes = len(modes_df) + modal_intensities = pd.DataFrame(index=range(n_modes)) + interacting_lasing_thresholds = np.inf * np.ones(n_modes) + candidates = [int(i) for i in np.where(lasing_thresholds < np.inf)[0]] + if not candidates: + return _finalise_modal_intensities( + modes_df, modal_intensities, interacting_lasing_thresholds + ) + + # Resolve the within-edge field for the hole burning: the bare-edge (per-edge + # mean) sampling over-clamps and spuriously suppresses co-lasing modes (it + # disagreed with Ge-Chong-Stone Eq. 28 on line_PRA, lasing one mode where two + # lase). ``oversample_size=None`` now auto-picks a wavelength-resolving size; + # pass 0 to force the old bare-edge behaviour. + if oversample_size is None: + oversample_size = _auto_oversample_size( + graph, modes_df, resolution=oversample_resolution, node_cap=oversample_node_cap + ) + work_graph = graph if not oversample_size else oversample_graph(graph, oversample_size) + pump = np.asarray(work_graph.graph["params"]["pump"], dtype=float) + pump_mask = _get_mask_matrices(work_graph.graph["params"])[1] + # small, fixed budget for each trust-region sub-solve (a local, warm-started + # solve converges in tens of evaluations; independent of params["max_steps"]) + max_steps = 30 + + # linear competition matrix only for the amplitude *unit* (diagonal) -- the + # active set itself is found self-consistently, not borrowed from it + t_linear = compute_mode_competition_matrix(work_graph, modes_df) + t_diag = np.array( + [abs(t_linear[i, i]) if abs(t_linear[i, i]) > 1e-12 else 1.0 for i in range(n_modes)] + ) + threshold_modes = modes_df["threshold_lasing_modes"].to_numpy() + + # NB: do not add per-mode "0 at its own threshold" baseline columns -- a + # mode's threshold is off the shared pump grid, so the *other* modes are then + # undefined (NaN -> 0) at that pump, putting a spurious dip in their curves. + # The grid already records every mode at every pump, with 0 below activation. + + mode_state: dict[int, np.ndarray] = {} + field_state: dict[int, np.ndarray] = {} + a_state: dict[int, float] = {} + unit_scale: dict[int, float] = {} + + # cap every (k, a) solve's k-window by the spacing of the *full* candidate + # set: a sparse active set on a dense spectrum must not wander across + # neighbouring roots (the per-set spacing alone cannot see them) + cand_ks = np.sort([float(from_complex(threshold_modes[c])[0]) for c in candidates]) + k_cap = float(np.clip(0.2 * np.diff(cand_ks).min(), 1e-6, 0.1)) if len(cand_ks) > 1 else None + + def _onset_amplitude(c, d0): + """Linear-slope amplitude estimate ``(D0 - thr)/(T_cc thr)`` for warm starts. + + Bootstraps and adds used to start at the 1e-3 floor, which sits in the + basin of the trivial ``a = 0`` root of the bounded solve (the same trap + the old measured onset probe dodged by probing at 1.2x threshold); + starting at the physical near-threshold estimate keeps the solve on the + lasing branch. + """ + thr_c = float(lasing_thresholds[c]) + slope = 1.0 / (t_diag[c] * thr_c) if (t_diag[c] > 0.0 and thr_c > 0.0) else 0.0 + return max(1e-3, slope * max(float(d0) - thr_c, 0.0)) + + def _init(i): + mode = np.asarray(from_complex(threshold_modes[i]), dtype=float) + field = _single_mode_field_intensity( + graph_with_pump(work_graph, float(lasing_thresholds[i])), mode, pump_mask + ) + mode_state[i], field_state[i], a_state[i] = mode, field, 0.0 + unit_scale[i] = _newton_onset_unit_scale( + work_graph, mode, field, float(lasing_thresholds[i]), t_diag[i] + ) + + def _solve_into(active_ids, d0, outer=6): + modes_out, fields_out, a_out, converged = _solve_active_set( + work_graph, + [mode_state[i] for i in active_ids], + [field_state[i] for i in active_ids], + [a_state[i] if a_state[i] > 1e-3 else _onset_amplitude(i, d0) for i in active_ids], + d0, + pump, + pump_mask, + max_steps, + seed, + outer=outer, + k_window_cap=k_cap, + ) + for j, i in enumerate(active_ids): + mode_state[i], field_state[i], a_state[i] = ( + modes_out[j], + fields_out[j], + float(a_out[j]), + ) + return converged + + def _advance_active_set(active_ids, d0_from, d0_to): + """Continue the active set ``d0_from -> d0_to`` along the physical branch. + + The frozen-field solve, fully relaxed in one coarse pump step, can walk + the background field off the physical (maximal-output) branch into a + spurious *lower-total* fixed point of the iteration -- far above + threshold near a mode crossing this reallocates amplitude between + competing modes and drops the summed intensity (the visible kink), even + though the true SALT total is monotone in pump. The cure is a true + continuation: short pump sub-steps with *light* field relaxation + (``outer=2``), so the field tracks the slowly-moving branch instead of + being free to converge to the distant spurious point. (A controlled + experiment confirmed ``outer=1`` and ``outer=2`` agree and give a smooth, + monotone crossing where full relaxation jumps discontinuously.) On a + well-separated spectrum the field barely moves and this is just the + normal solve at finer resolution -- the Ge-Chong-Stone line_PRA result + is unchanged. + """ + n_sub = max(1, int(np.ceil(abs(d0_to - d0_from) / 0.02))) + for d0 in np.linspace(d0_from, d0_to, n_sub + 1)[1:]: + _solve_into(active_ids, float(d0), outer=2) + + active: list[int] = [] # confirmed lasing ids, carried along the continuation + first = float(np.min(lasing_thresholds[candidates])) + d0_prev = None + prev_state: dict[int, tuple] = {} # last accepted (mode, field, a) per active id + prev_total = 0.0 # last accepted summed (unit-scaled) output -- a monotone floor + for D0 in np.linspace(first, max_pump_intensity, D0_steps): + if active and d0_prev is not None: + # pump-step continuity: advance the carried-over set in-basin + _advance_active_set(list(active), d0_prev, float(D0)) + d0_prev = float(D0) + bootstrapped_off: set[int] = set() # bootstrapped then vanished at this D0 + rejected_adds: set[int] = set() # adds reverted by the continuity guard at this D0 + pending_add = None # candidate added on the previous sweep, not yet vetted + veteran_snapshot: dict[int, tuple] = {} + pre_add_total = 0.0 # summed (unit-scaled) output before the pending add + for _ in range(3 * len(candidates) + 2): # active-set sweeps until stable + if active: + modes_out, fields_out, a_out, converged = _solve_active_set( + work_graph, + [mode_state[i] for i in active], + [field_state[i] for i in active], + [a_state[i] if a_state[i] > 1e-3 else _onset_amplitude(i, D0) for i in active], + D0, + pump, + pump_mask, + max_steps, + seed, + k_window_cap=k_cap, + ) + if not converged: + warnings.warn( + f"full_salt_newton field loop did not fully converge at D0={D0:.4g}.", + stacklevel=2, + ) + for j, i in enumerate(active): + mode_state[i], field_state[i], a_state[i] = ( + modes_out[j], + fields_out[j], + float(a_out[j]), + ) + if pending_add is not None: + # Continuity guard: SALT total output is monotone in the pump, + # so an *add* whose re-solve *lowers* the summed intensity is a + # wrong-basin solution -- physically a newcomer turning on can + # only raise the total (it adds output and steals at most a + # little from the others, which nearly cancels: Ge-Chong-Stone). + # On a dense spectrum the coupled solve can instead flip to a + # different fixed point where the newcomer dominates and a + # stronger veteran is suppressed -- a partial collapse that + # drops the total and reads as a spurious kink (observed on the + # mini-buffon example, dk = 0.004). Total-monotonicity catches + # those partial collapses that a per-veteran "killed outright" + # test misses. Revert the veterans, reject the newcomer at this + # pump step (retried at the next), and re-probe. + post_total = sum( + a_state[i] * unit_scale.get(i, 1.0) for i in active if a_state[i] > 1e-4 + ) + # (a) total-dropping partial collapse, and (b) a total- + # *preserving* identity swap where a higher-threshold newcomer + # kills a lower-threshold veteran outright (the twin takeover: + # two modes closer than the solve can resolve in amplitude, so + # the indeterminate split dumps all of it on the newcomer -- + # the total is unchanged, so (a) alone misses it). Either is a + # wrong basin; keep the cluster on the veteran we were tracking. + killed_lower = any( + a_prev > 1e-2 + and a_state[i] <= 1e-3 + and lasing_thresholds[i] < lasing_thresholds[pending_add] + for i, (_m, _f, a_prev) in veteran_snapshot.items() + ) + wrong_basin = killed_lower or post_total < (1.0 - 1e-2) * pre_add_total + if wrong_basin: + for i, (m_prev, f_prev, a_prev) in veteran_snapshot.items(): + mode_state[i], field_state[i], a_state[i] = m_prev, f_prev, a_prev + active.remove(pending_add) + a_state[pending_add] = 0.0 + rejected_adds.add(pending_add) + pending_add = None + veteran_snapshot = {} + continue + if a_state[pending_add] <= 1e-4: + rejected_adds.add(pending_add) # stillborn add: do not retry at this D0 + pending_add = None + veteran_snapshot = {} + dropped = [i for i in active if a_state[i] <= 1e-4] + bootstrapped_off.update(dropped) + active = [i for i in active if a_state[i] > 1e-4] # drop vanished + if not active: + # Bootstrap an empty set: on the unsaturated background the + # noninteracting threshold is exact, and the gain-crossing probe + # below would wrongly reject a mode sitting *at* its threshold + # (alpha = 0 there, not < 0 -- so the first mode would otherwise + # turn on a full pump step late). Activate the lowest-threshold + # eligible candidate directly. + eligible = [ + c + for c in candidates + if lasing_thresholds[c] <= D0 and c not in bootstrapped_off + ] + if not eligible: + break + c = min(eligible, key=lambda i: lasing_thresholds[i]) + if c not in mode_state: + _init(c) + a_state[c] = _onset_amplitude(c, D0) + active.append(c) + continue # solve the bootstrapped set before probing others + # gain the not-yet-lasing candidates see on the current saturated background + background = _saturated_graph_multi( + work_graph, + [mode_state[i] for i in active], + [a_state[i] for i in active], + D0, + pump, + [field_state[i] for i in active], + ) + active_ks = [float(mode_state[i][0]) for i in active] + best, best_alpha, best_k = None, -1e-6, 0.0 + for c in candidates: + if c in active or c in rejected_adds or lasing_thresholds[c] > D0: + continue + if c not in mode_state: + _init(c) + # probe window consistent with the solve's spacing-tracking + # k-window (a fixed wide window lets the probe wander to a + # different branch than the solve will then confine it to) + gap = min(abs(float(mode_state[c][0]) - k) for k in active_ks) + window = float(np.clip(0.2 * gap, 1e-6, 0.3)) + kc = _refine_local( + mode_state[c], background, 1e-9, max_steps, seed, k_window=window + ) + # alpha = mode[1] < 0 => net gain => the mode crossed its + # *interacting* threshold. Tight margin (-1e-6): gain clamping + # pins an above-threshold mode's alpha at ~0^-, so a looser + # cutoff adds the mode pump-steps late and snaps its L--I curve; + # the a < 1e-4 drop rule above is the safety net against false + # adds. Add only the *most* above-threshold candidate per sweep: + # several at once hand the coupled solve a multistable warm + # start, which can converge to the wrong basin and zero the true + # winner (observed on line_PRA, where a simultaneous three-mode + # add transiently deleted the dominant mode). + if kc[1] < best_alpha: + best, best_alpha, best_k = c, float(kc[1]), float(kc[0]) + if best is None: + break + veteran_snapshot = { + i: (mode_state[i].copy(), np.asarray(field_state[i]).copy(), a_state[i]) + for i in active + } + pre_add_total = sum(a_state[i] * unit_scale.get(i, 1.0) for i in active) + pending_add = best + mode_state[best] = np.array([best_k, 0.0]) + # A newcomer crosses its *interacting* threshold here, so it turns on + # from ~0 -- warm-start it small. The bare-threshold onset estimate + # (right for the empty-set bootstrap, where there is no suppression) + # over-shoots badly when D0 >> bare threshold and pulls the coupled + # solve into a wrong basin where the newcomer dominates a stronger + # veteran. Cap it well below the established amplitudes; the field is + # anchored by the veterans, so the trivial a = 0 root is not a risk. + a_cap = 1e-2 * max(a_state[i] for i in active) + a_state[best] = float(min(_onset_amplitude(best, D0), max(a_cap, 1e-3))) + active.append(best) + # Total-output ratchet (physical floor). A steady-state laser's total + # output cannot fall as the pump rises. Far above threshold (~>15x) the + # frozen-field single-pole iteration can converge to a spurious + # lower-total branch -- a mode reallocation it renders as an unphysical + # dip -- robustly across relaxation / step size (a limit of the method, + # not a tuning bug). Enforce monotonicity directly: if the total dropped, + # hold any veteran whose amplitude collapsed back to its last accepted + # value (new modes still join and raise the total). The curve then + # plateaus where this bites -- an honest marker of the per-mode + # resolution limit -- instead of kinking downward. + cur_total = sum(a_state.get(i, 0.0) * unit_scale.get(i, 1.0) for i in active) + if prev_state and cur_total < (1.0 - 1e-2) * prev_total: + for i in list(active): + if i in prev_state and a_state[i] < 0.5 * prev_state[i][2]: + m_prev, f_prev, a_prev = prev_state[i] + mode_state[i], field_state[i], a_state[i] = m_prev.copy(), f_prev.copy(), a_prev + warnings.warn( + f"full_salt_newton: held collapsing mode(s) at D0={D0:.4g} to keep the total " + "output monotone (per-mode resolution limit far above threshold).", + stacklevel=2, + ) + prev_total = max( + prev_total, sum(a_state.get(i, 0.0) * unit_scale.get(i, 1.0) for i in active) + ) + prev_state = { + i: (mode_state[i].copy(), np.asarray(field_state[i]).copy(), a_state[i]) for i in active + } + + for i in candidates: + value = max(a_state.get(i, 0.0) * unit_scale.get(i, 1.0), 0.0) if i in active else 0.0 + modal_intensities.loc[i, D0] = value + if i in active and a_state[i] > 0 and D0 < interacting_lasing_thresholds[i]: + interacting_lasing_thresholds[i] = D0 + + return _finalise_modal_intensities(modes_df, modal_intensities, interacting_lasing_thresholds) + + def pump_trajectories(modes_df, graph, return_approx=False, quality_method="eigenvalue"): """For a sequence of D0s, find the mode positions of the modes modes.""" diff --git a/netsalt/params.py b/netsalt/params.py index 1e20ac0f..2ef004ac 100644 --- a/netsalt/params.py +++ b/netsalt/params.py @@ -100,6 +100,33 @@ class NetSaltParams(BaseModel): reduction_factor: float | None = None n_modes_max: int | None = None + # --- Modal-intensity solver -------------------------------------------- + # ``intensity_method`` selects how the lasing L--I curves are computed from + # the threshold modes and the mode-competition matrix. Accepted values: + # ``"linear"`` (default) — the near-threshold SALT model: one + # pump-independent competition matrix and a linear solve + # (``compute_modal_intensities``). Piecewise-linear curves; fast. + # ``"full_salt_newton"`` — operator-level SALT: solves the saturated + # nonlinear eigenproblem for the active modes' ``(k, a)`` (real k, + # amplitude), capturing gain-clamping mode suppression. Validated against + # Ge-Chong-Stone on examples/line_PRA; slower. + # See doc/source/lasing.rst and issue #42. + intensity_method: Literal["linear", "full_salt_newton"] | None = None + # full_salt_newton knobs: + intensity_max_iter: int | None = None + intensity_tol: float | None = None + intensity_damping: float | None = None + salt_D0_steps: int | None = None + # edge_size passed to oversample_graph for the within-edge hole-burning + # resolution (None auto-picks a wavelength-resolving size). + intensity_oversample_size: float | None = None + # full_salt_newton auto-oversampling knobs: sub-edges per wavelength + # (resolution) and the node cap that bounds the oversampled operator. The cap + # makes the solver feasible on large graphs (buffon) but reduces accuracy + # there; raise it (slower) when the modal magnitudes matter. + intensity_oversample_resolution: int | None = None + intensity_oversample_node_cap: int | None = None + # --- Infrastructure ---------------------------------------------------- n_workers: int | None = None diff --git a/netsalt/physics.py b/netsalt/physics.py index 0580ed4a..6fc6ac21 100644 --- a/netsalt/physics.py +++ b/netsalt/physics.py @@ -137,6 +137,39 @@ def dispersion_relation_pump(freq, params=None): ) +def dispersion_relation_pump_saturated(freq, params=None): + r"""Saturated pumped dispersion relation (full SALT, issue #42). + + Identical to :func:`dispersion_relation_pump` except the scalar gain term + :math:`D_0\,\delta_\mathrm{pump}` is replaced by a per-edge *saturated* + effective-pump field ``params["D0_eff"]`` carrying the spatial-hole-burning + denominator :math:`1 + \sum_\nu \Gamma_\nu a_\nu |\Psi_\nu|^2`: + + .. math:: + + k(\omega) = \frac{\omega}{c} \sqrt{\epsilon + \gamma(\omega)\, D_0^\mathrm{eff}}. + + With ``D0_eff = D0 * pump`` (the unsaturated limit, denominator one) this + reduces *exactly* to :func:`dispersion_relation_pump`, so it can be swapped + in without changing passive behaviour. + + Args: + freq (float): frequency + params (dict): parameters, must include ``dielectric_constant``; if + ``D0_eff`` (a per-edge array) is absent it falls back to + :func:`dispersion_relation_pump`. + """ + if not params: + raise ValueError("Please provide dispersion parameters") + + dielectric = np.asarray(params["dielectric_constant"]) + c = params.get("c", 1.0) + if "D0_eff" not in params: + return dispersion_relation_pump(freq, params) + + return freq * np.sqrt(dielectric + gamma(freq, params) * np.asarray(params["D0_eff"])) / c + + def set_dielectric_constant(graph, params, custom_values=None, rng=None): """Set dielectric constant in params, from dielectric constant or refraction index. diff --git a/netsalt/pipeline.py b/netsalt/pipeline.py index 8c87d6dc..8870ac42 100644 --- a/netsalt/pipeline.py +++ b/netsalt/pipeline.py @@ -42,6 +42,7 @@ ) from .modes import ( compute_modal_intensities, + compute_modal_intensities_full_salt_newton, compute_mode_competition_matrix, find_passive_modes, find_threshold_lasing_modes, @@ -354,9 +355,14 @@ def step_compute_mode_competition_matrix( def step_compute_modal_intensities( - p: NetSaltParams, threshold_modes_df, competition_matrix, lasing_modes_id + p: NetSaltParams, qg, threshold_modes_df, competition_matrix, pump, lasing_modes_id ): - """Compute modal intensities over the pump-strength sweep.""" + """Compute modal intensities over the pump-strength sweep. + + Dispatches on ``params["intensity_method"]`` (default ``"linear"``). The + ``full_salt_newton`` solver needs the pumped graph, so it is reattached + here; the ``linear`` path is unchanged and ignores it. + """ out = _outdir(p) / _apply_lasing_ids("modal_intensities.h5", lasing_modes_id) if out.exists() and not _force(p): return load_modes(str(out)) @@ -365,7 +371,30 @@ def step_compute_modal_intensities( D0_max = p.get("intensities_D0_max") if D0_max is None: D0_max = p.get("D0_max", 0.1) - modes_df = compute_modal_intensities(threshold_modes_df, D0_max, competition_matrix) + + method = p.get("intensity_method") or "linear" + if method != "linear": + # the pump-dependent solvers evaluate profiles on the pumped graph + qg = _attach_pump_to_graph(p, qg, pump) + if method == "linear": + modes_df = compute_modal_intensities(threshold_modes_df, D0_max, competition_matrix) + elif method == "full_salt_newton": + modes_df = compute_modal_intensities_full_salt_newton( + qg, + threshold_modes_df, + D0_max, + D0_steps=p.get("salt_D0_steps", 30), + tol=p.get("intensity_tol", 1e-8), + oversample_size=p.get("intensity_oversample_size"), + inner_max_iter=p.get("intensity_max_iter", 25), + inner_damping=p.get("intensity_damping", 0.8), + oversample_resolution=p.get("intensity_oversample_resolution", 12), + oversample_node_cap=p.get("intensity_oversample_node_cap", 3000), + ) + else: # pragma: no cover - guarded by the NetSaltParams Literal + raise ValueError( + f"Unknown intensity_method {method!r}; expected 'linear' or 'full_salt_newton'." + ) save_modes(modes_df, filename=str(out)) return modes_df @@ -650,7 +679,7 @@ def compute_lasing_modes(p: NetSaltParams, lasing_modes_id=None): plot_mode_competition_matrix_fig(p, competition, lasing_modes_id) intensities_df = step_compute_modal_intensities( - p, threshold_modes_df, competition, lasing_modes_id + p, qg, threshold_modes_df, competition, pump, lasing_modes_id ) plot_ll_curve_fig(p, qg, intensities_df, lasing_modes_id) plot_stem_spectra_fig(p, qg, intensities_df, lasing_modes_id) @@ -684,7 +713,9 @@ def compute_controllability(p: NetSaltParams): trajectories_df = step_compute_mode_trajectories(p, qg, passive_modes_df, pump, ids) threshold_modes_df = step_find_threshold_modes(p, qg, trajectories_df, pump, ids) competition = step_compute_mode_competition_matrix(p, qg, threshold_modes_df, pump, ids) - intensities_df = step_compute_modal_intensities(p, threshold_modes_df, competition, ids) + intensities_df = step_compute_modal_intensities( + p, qg, threshold_modes_df, competition, pump, ids + ) plot_ll_curve_fig(p, qg, intensities_df, ids) spectra = intensities_df[ diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index 26c215fb..2e3ed0bb 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -21,6 +21,21 @@ L = logging.getLogger(__name__) +# Quantum-graph laplacians at or below this dimension use a direct dense +# eigensolve instead of ARPACK shift-invert: ``eigs(sigma=0)`` carries a fixed +# per-call overhead (a sparse LU factorisation + Arnoldi restart, ~2.5 ms) that +# only dominates for small graphs, where ``np.linalg.eig`` on the dense matrix is +# faster and returns the same nearest-zero eigenpair. Above the crossover (~50 +# nodes, measured) ARPACK wins decisively: it is ~flat in N on the banded +# quantum-graph laplacian (~2.5 ms at N=60, ~6 ms at N=1000), whereas dense is +# O(N^3) (~24 ms at N=120, ~160 ms at N=250). Kept high (256) here for robust +# *mode finding* on dense-spectrum 2D graphs (e.g. buffon): the grid scan probes +# many near-singular k, where ARPACK shift-invert (sigma=0, an LU of a near-singular +# matrix) is slow/unstable, while dense is robust. ``full_salt_newton`` lowers it +# *locally* (NEWTON_DENSE_EIG_MAX) for its own banded, oversampled saturated solves, +# where ARPACK is both fast (~flat in N) and stable (isolated lasing modes). +DENSE_EIG_MAX = 256 + def create_quantum_graph( graph, params=None, positions=None, lengths=None, seed=42, noise_level=0.001 @@ -507,6 +522,21 @@ def laplacian_quality(laplacian, method="eigenvalue", rng=None): starting vector. If None, a fresh generator with fresh entropy is created. Pass a seeded generator for reproducibility. """ + # Dense fast path for small matrices (see DENSE_EIG_MAX): the nearest-zero + # eigenvalue is the smallest-magnitude one, identical to ``eigs(sigma=0)`` but + # without ARPACK's per-call overhead. ``rng`` is irrelevant here (no ARPACK + # start vector), keeping the result deterministic. + if method in ("eigenvalue", "complex_eigenvalue") and laplacian.shape[0] <= DENSE_EIG_MAX: + dense = laplacian.toarray() + # a root-finder can probe a ``k`` whose operator overflows to inf/NaN; + # ``np.linalg.eigvals`` raises there, so signal "not a mode" (quality 1) + # exactly as the ARPACK branch does on non-convergence. + if not np.isfinite(dense).all(): + return 1.0 if method == "eigenvalue" else 1.0 + 0j + eigenvalues = np.linalg.eigvals(dense) + lam = eigenvalues[np.argmin(np.abs(eigenvalues))] + return abs(lam) if method == "eigenvalue" else complex(lam) + if rng is None: rng = np.random.default_rng() v0 = rng.random(laplacian.shape[0]) diff --git a/tests/test_functional.py b/tests/test_functional.py index bca019fd..db84e187 100644 --- a/tests/test_functional.py +++ b/tests/test_functional.py @@ -87,3 +87,183 @@ def test_ComputeLasingModes(working_directory): assert_equal_trees( expected_dir, result_dir, specific_args={"out": {"patterns": [r".*\.h5$"], "atol": 1e-5}} ) + + +def test_ComputeLasingModes_full_salt_newton(working_directory): + """End-to-end smoke test of the operator-level Newton solver. + + Runs the same pipeline with ``intensity_method="full_salt_newton"`` on a + coarse pump grid: it exercises the multimode driver, the amplitude solve and + the frequency/profile fixed point on genuine threshold modes (no byte + reference -- this is a coverage/regression smoke test). The dominant mode must + lase below the max pump, matching the linear model's lowest threshold. + """ + create_graph() + + params = load_config("config.yaml") + params["intensity_method"] = "full_salt_newton" + params["salt_D0_steps"] = 2 + compute_lasing_modes(params) + + result_dir, _ = working_directory + intensities = netsalt.load_modes(str(result_dir / "modal_intensities.h5")) + intensity_cols = [c for c in intensities.columns if c[0] == "modal_intensities"] + assert intensity_cols, "no modal-intensity columns were written" + # at least one mode reaches a positive intensity at the largest pump + pumps = sorted(c[1] for c in intensity_cols) + at_max = intensities[("modal_intensities", pumps[-1])].to_numpy(dtype=float) + assert np.nan_to_num(at_max).max() > 0.0 + + +def test_full_salt_newton_multimode_two_ring(): + """full_salt_newton lases *several* modes on a genuinely multimode graph. + + Two detuned rings (different sizes) joined by a bridge: the detuning localises + each mode onto one ring, so with a narrow gain they barely compete and several + co-lase. Guards the multimode path -- in particular against the collapse where + a loose k-window let the solve drop a co-lasing mode by drifting its k to a + spurious a=0 root (it then reported a single, wrong mode). + """ + import networkx as nx + + from netsalt.modes import ( + compute_modal_intensities_full_salt_newton, + find_passive_modes, + find_threshold_lasing_modes, + pump_trajectories, + ) + from netsalt.physics import dispersion_relation_pump + from netsalt.quantum_graph import create_quantum_graph, set_total_length + + n_a, n_b = 7, 9 + g = nx.disjoint_union(nx.cycle_graph(n_a), nx.cycle_graph(n_b)) + g.add_edge(0, n_a) + g.add_edge(2, n_a + n_b) + g.add_edge(n_a + 4, n_a + n_b + 1) + pos = {} + for i in range(n_a): + pos[i] = [-1.6 + 0.9 * np.cos(2 * np.pi * i / n_a), 0.9 * np.sin(2 * np.pi * i / n_a)] + for j in range(n_b): + pos[n_a + j] = [ + 1.7 + 1.25 * np.cos(2 * np.pi * j / n_b), + 1.25 * np.sin(2 * np.pi * j / n_b), + ] + pos[n_a + n_b] = [-1.6, -2.2] + pos[n_a + n_b + 1] = [1.7, -2.6] + positions = np.array([pos[i] for i in range(len(g))]) + params = { + "open_model": "open", + "c": 1.0, + "k_a": 3.567, + "gamma_perp": 0.12, + "k_min": 3.45, + "k_max": 3.66, + "alpha_min": -0.05, + "alpha_max": 0.15, + "n_workers": 1, + "n_modes_max": 10, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 1.0, + "D0_steps": 14, + "dielectric_params": { + "method": "uniform", + "inner_value": 9.0, + "outer_value": 1.0, + "loss": 0.0, + }, + } + create_quantum_graph(g, params, positions=positions) + set_total_length(g, 9.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + + passive = find_passive_modes(g, method="contour") + assert len(passive) >= 2, "fixture should find several modes" + g.graph["params"]["pump"] = np.array([1.0 if g[u][v]["inner"] else 0.0 for u, v in g.edges()]) + trajectories = pump_trajectories(passive, g, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, g) + + df = compute_modal_intensities_full_salt_newton(g, tdf.copy(), 1.0, D0_steps=5) + cols = [c for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + data = np.nan_to_num(df[cols].to_numpy(dtype=float)) + peak = max(data.max(), 1e-9) + n_lasing = int(np.sum(data.max(axis=1) > 1e-2 * peak)) + assert n_lasing >= 2, f"expected multimode lasing, got {n_lasing}" + + +def test_full_salt_newton_multimode_chaotic_ring(): + """full_salt_newton lases several modes on a single ring with random chords. + + A 14-node ring plus six fixed chords (the buffon mechanism shrunk down): the + chords close extra loops, so the spectrum is dense and the modes localise on + different loops. With a narrow gain on a mode cluster, several co-lase on one + small single-component graph -- guards the multimode example + ``examples/chaotic_ring``. + """ + import networkx as nx + + from netsalt.modes import ( + compute_modal_intensities_full_salt_newton, + find_passive_modes, + find_threshold_lasing_modes, + pump_trajectories, + ) + from netsalt.physics import dispersion_relation_pump + from netsalt.quantum_graph import create_quantum_graph, set_total_length + + m = 14 + chords = [(0, 7), (1, 5), (2, 12), (3, 5), (3, 10), (7, 11)] + g = nx.cycle_graph(m) + g.add_edges_from(chords) + g.add_edge(0, m) + g.add_edge(m // 2, m + 1) + pos = {i: [np.cos(2 * np.pi * i / m), np.sin(2 * np.pi * i / m)] for i in range(m)} + pos[m] = [1.6, 0.0] + pos[m + 1] = [-1.6, 0.0] + positions = np.array([pos[i] for i in range(len(g))]) + params = { + "open_model": "open", + "c": 1.0, + "k_a": 2.81, + "gamma_perp": 0.10, + "k_min": 2.55, + "k_max": 3.25, + "alpha_min": -0.05, + "alpha_max": 0.25, + "n_workers": 1, + "n_modes_max": 40, + "quality_threshold": 1e-3, + "search_stepsize": 0.005, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "D0_max": 0.5, + "D0_steps": 14, + "dielectric_params": { + "method": "uniform", + "inner_value": 9.0, + "outer_value": 1.0, + "loss": 0.0, + }, + } + create_quantum_graph(g, params, positions=positions) + set_total_length(g, 12.0) + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + + passive = find_passive_modes(g, method="contour") + assert len(passive) >= 3, "fixture should find several modes" + g.graph["params"]["pump"] = np.array([1.0 if g[u][v]["inner"] else 0.0 for u, v in g.edges()]) + trajectories = pump_trajectories(passive, g, return_approx=True) + tdf = find_threshold_lasing_modes(trajectories, g) + + df = compute_modal_intensities_full_salt_newton(g, tdf.copy(), 0.5, D0_steps=5) + cols = [c for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + data = np.nan_to_num(df[cols].to_numpy(dtype=float)) + peak = max(data.max(), 1e-9) + n_lasing = int(np.sum(data.max(axis=1) > 1e-2 * peak)) + assert n_lasing >= 3, f"expected multimode lasing, got {n_lasing}" diff --git a/tests/test_unit.py b/tests/test_unit.py index a238450e..f3ed601a 100644 --- a/tests/test_unit.py +++ b/tests/test_unit.py @@ -1501,3 +1501,392 @@ def test_threshold_in_middle_column(self): # |imag| minimal in the middle column -> +1 stays in range. df = self._modes_df([1.0, 0.0, 1.0]) plot_pump_traj(df) + + +class TestSaturatedDispersion: + """``dispersion_relation_pump_saturated`` (full-SALT gain term).""" + + def _params(self): + return { + "dielectric_constant": np.array([2.0, 3.0, 2.5]), + "pump": np.array([1.0, 0.0, 1.0]), + "D0": 0.03, + "c": 1.0, + "gamma_perp": 0.5, + "k_a": 5.0, + } + + def test_reduces_to_pumped_when_unsaturated(self): + """With ``D0_eff = D0 * pump`` (denominator one) it must reproduce + ``dispersion_relation_pump`` exactly.""" + from netsalt.physics import ( + dispersion_relation_pump, + dispersion_relation_pump_saturated, + ) + + params = self._params() + base = dispersion_relation_pump(5.1, params) + sat_params = dict(params, D0_eff=params["D0"] * params["pump"]) + assert np.allclose(base, dispersion_relation_pump_saturated(5.1, sat_params)) + + def test_falls_back_to_pumped_without_D0_eff(self): + from netsalt.physics import ( + dispersion_relation_pump, + dispersion_relation_pump_saturated, + ) + + params = self._params() + assert np.allclose( + dispersion_relation_pump(5.1, params), + dispersion_relation_pump_saturated(5.1, params), + ) + + def test_saturation_lowers_the_gain_contribution(self): + """A larger denominator (D0_eff < D0*pump) pulls k toward the passive + dielectric value on the pumped edges.""" + from netsalt.physics import ( + dispersion_relation_dielectric, + dispersion_relation_pump, + dispersion_relation_pump_saturated, + ) + + params = self._params() + passive = dispersion_relation_dielectric(5.1, params) + pumped = dispersion_relation_pump(5.1, params) + sat_params = dict(params, D0_eff=0.5 * params["D0"] * params["pump"]) + saturated = dispersion_relation_pump_saturated(5.1, sat_params) + # on the pumped edges the saturated k sits between passive and full pump + gain_edges = params["pump"] > 0 + assert np.all( + np.abs(saturated - passive)[gain_edges] < np.abs(pumped - passive)[gain_edges] + ) + + +class TestModeOnNodesQualityFlag: + """``check_quality=False`` lets the profile helpers evaluate off-threshold.""" + + def _line_graph(self, **kw): + return make_line_graph(**kw) + + def test_check_quality_false_returns_vector_on_non_mode(self): + from netsalt.modes import mode_on_nodes + + g = self._line_graph(n_edges=4) + g.graph["params"]["quality_threshold"] = 1e-12 + # would raise with the default check; must not with it disabled + sol = mode_on_nodes([3.0, 0.05], g, check_quality=False) + assert sol.shape == (len(g),) + + def test_mean_mode_on_edges_threads_the_flag(self): + import netsalt + from netsalt.modes import mean_mode_on_edges + from netsalt.physics import dispersion_relation_pump + from netsalt.quantum_graph import update_parameters + + g = self._line_graph(n_edges=4) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + update_parameters( + g, {"k_a": 3.0, "gamma_perp": 1.0, "D0": 0.5, "pump": np.ones(len(g.edges))} + ) + g.graph["params"]["quality_threshold"] = 1e-12 + mean = mean_mode_on_edges([3.0, 0.05], g, check_quality=False) + assert mean.shape == (len(g.edges),) + + +class TestIntensitySolveHelpers: + """Pure-algebra helpers shared by the intensity solvers.""" + + def test_slopes_shifts_identity_matrix(self): + from netsalt.modes import _intensity_slopes_shifts + + thresholds = np.array([2.0, 4.0]) + T = np.eye(2) + slopes, shifts = _intensity_slopes_shifts(T, thresholds, [0, 1]) + # T = I -> slopes = 1/threshold, shifts = 1, so intensity(D0) = D0/thr - 1 + assert np.allclose(slopes, 1.0 / thresholds) + assert np.allclose(shifts, 1.0) + + def test_finalise_writes_sorted_intensity_columns(self): + import pandas as pd + + from netsalt.modes import _finalise_modal_intensities + + modal = pd.DataFrame(index=range(2)) + modal.loc[0, 0.5] = 0.0 + modal.loc[0, 0.2] = 0.0 # inserted out of order on purpose + modal.loc[0, 0.8] = 1.0 + out = _finalise_modal_intensities( + pd.DataFrame(index=range(2)), modal, np.array([0.5, np.inf]) + ) + pumps = [c[1] for c in out.columns if c[0] == "modal_intensities"] + assert pumps == sorted(pumps) + assert np.allclose(out["interacting_lasing_thresholds"].to_numpy(), [0.5, np.inf]) + + +class TestIntensityMethodDispatch: + """``step_compute_modal_intensities`` routes on ``intensity_method``.""" + + def _params(self, tmp_path, method): + from netsalt.params import NetSaltParams + + extra = {} if method is None else {"intensity_method": method} + return NetSaltParams.from_dict( + {"outdir": str(tmp_path), "force": True, "intensities_D0_max": 1.0, **extra} + ) + + def _run(self, tmp_path, monkeypatch, method): + import pandas as pd + + from netsalt import pipeline + + calls = [] + + def make(name): + def _fake(*args, **kwargs): + calls.append(name) + return pd.DataFrame({("modal_intensities", 0.5): [0.0]}) + + return _fake + + monkeypatch.setattr(pipeline, "compute_modal_intensities", make("linear")) + monkeypatch.setattr( + pipeline, "compute_modal_intensities_full_salt_newton", make("full_salt_newton") + ) + monkeypatch.setattr(pipeline, "_attach_pump_to_graph", lambda p, qg, pump: qg) + monkeypatch.setattr(pipeline, "save_modes", lambda *a, **k: None) + + p = self._params(tmp_path, method) + pipeline.step_compute_modal_intensities( + p, object(), pd.DataFrame(), np.zeros((1, 1)), None, None + ) + return calls + + def test_default_is_linear(self, tmp_path, monkeypatch): + assert self._run(tmp_path, monkeypatch, None) == ["linear"] + + def test_dispatches_each_method(self, tmp_path, monkeypatch): + for method in ("linear", "full_salt_newton"): + assert self._run(tmp_path, monkeypatch, method) == [method] + + +def _independent_lasing_fixture(): + """Build a small open dielectric line cavity and return ``(graph, threshold_df)``. + + A short Fabry--Perot straddling the gain line at ``k_a = 15`` with a handful of + well-separated lasing modes -- an independent fixture (not ``line_PRA``) for + the full-SALT consistency tests. Runs the real passive -> pump -> trajectories + -> threshold pipeline so the modes are genuine. + """ + import networkx as nx + + import netsalt + from netsalt.modes import find_threshold_lasing_modes, pump_trajectories, scan_frequencies + from netsalt.physics import dispersion_relation_pump + from netsalt.quantum_graph import create_quantum_graph, set_total_length + + n_edges = 8 + g = nx.path_graph(n_edges + 1) + positions = np.array([[float(i), 0.0] for i in range(n_edges + 1)]) + params = { + "open_model": "open", + "c": 1.0, + "k_a": 15.0, + "gamma_perp": 3.0, + "k_min": 12.0, + "k_max": 18.0, + "k_n": 80, + "alpha_min": 0.0, + "alpha_max": 1.0, + "alpha_n": 20, + "quality_threshold": 1e-3, + "search_stepsize": 0.01, + "max_steps": 1000, + "max_tries_reduction": 50, + "reduction_factor": 0.8, + "n_workers": 1, + "D0_max": 1.0, + "D0_steps": 10, + "dielectric_params": { + "method": "uniform", + "inner_value": 9.0, + "outer_value": 1.0, + "loss": 0.0, + }, + } + create_quantum_graph(g, params, positions=positions) + set_total_length(g, 0.5) # short -> well-separated longitudinal modes + netsalt.set_dielectric_constant(g, g.graph["params"]) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + + qualities = scan_frequencies(g) + passive = netsalt.find_passive_modes( + g, qualities, method="grid", min_distance=2, threshold_abs=0.1 + ) + pump = np.array([1.0 if g[u][v]["inner"] else 0.0 for u, v in g.edges()]) + g.graph["params"]["pump"] = pump + trajectories = pump_trajectories(passive, g, return_approx=True) + return g, find_threshold_lasing_modes(trajectories, g) + + +class TestFullSaltNewton: + """Building blocks of the operator-level single-mode Newton solver.""" + + def _pump_graph(self, n_edges=4): + import netsalt + from netsalt.physics import dispersion_relation_pump + from netsalt.quantum_graph import update_parameters + + g = make_line_graph(n_edges=n_edges) + netsalt.set_dispersion_relation(g, dispersion_relation_pump) + update_parameters( + g, {"k_a": 3.0, "gamma_perp": 1.0, "D0": 0.0, "pump": np.ones(len(g.edges))} + ) + g.graph["params"]["quality_threshold"] = 10.0 + return g + + def test_intensity_method_literal_accepts_newton(self): + from netsalt.params import NetSaltParams + + assert NetSaltParams(intensity_method="full_salt_newton")["intensity_method"] == ( + "full_salt_newton" + ) + + def test_saturated_graph_reduces_to_pumped_at_zero_amplitude(self): + from netsalt.modes import _saturated_graph_at + from netsalt.physics import dispersion_relation_pump_saturated + + g = self._pump_graph() + pump = np.asarray(g.graph["params"]["pump"], dtype=float) + field = np.ones(len(g.edges)) + gsat = _saturated_graph_at(g, [3.0, 0.0], 0.0, 0.5, pump, field) + # a = 0 -> denominator 1 -> D0_eff = D0 * pump (unsaturated), saturated + # dispersion swapped in. (Equivalence to dispersion_relation_pump at this + # D0_eff is covered by TestSaturatedDispersion.) + np.testing.assert_allclose(gsat.graph["params"]["D0_eff"], 0.5 * pump) + assert gsat.graph["dispersion_relation"] is dispersion_relation_pump_saturated + # the throwaway copy must not mutate the original graph + assert "D0_eff" not in g.graph["params"] + + def test_saturated_graph_lowers_effective_pump_with_amplitude(self): + from netsalt.modes import _saturated_graph_at + + g = self._pump_graph() + pump = np.asarray(g.graph["params"]["pump"], dtype=float) + field = np.ones(len(g.edges)) + unsat = _saturated_graph_at(g, [3.0, 0.0], 0.0, 0.5, pump, field).graph["params"]["D0_eff"] + sat = _saturated_graph_at(g, [3.0, 0.0], 1.0, 0.5, pump, field).graph["params"]["D0_eff"] + # hole burning reduces the effective pump on the gain edges + assert np.all(sat[pump > 0] < unsat[pump > 0]) + + def test_field_intensity_is_finite_per_edge(self): + from netsalt.modes import _get_mask_matrices, _single_mode_field_intensity + + g = self._pump_graph() + pump_mask = _get_mask_matrices(g.graph["params"])[1] + inten = _single_mode_field_intensity(g, [3.0, 0.05], pump_mask) + assert inten.shape == (len(g.edges),) + assert np.all(np.isfinite(inten)) + + def test_reduces_to_linear_onset_slope_on_independent_graph(self): + """full_salt_newton's reported intensity is in the linear modal-intensity + unit on a graph *other* than line_PRA: the dominant mode's onset slope + matches the linear ``1/(T_μμ·D0_thr)``. Guards the unit-consistency fix + against the graph-dependent within-edge form factor.""" + from netsalt.modes import ( + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + ) + + g, tdf = _independent_lasing_fixture() + + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + assert np.any(thresholds < np.inf), "fixture must produce a lasing mode" + t0 = int(np.argmin(thresholds)) + thr0 = float(thresholds[t0]) + T = compute_mode_competition_matrix(g, tdf) + linear_slope = 1.0 / (T[t0, t0] * thr0) + + # measure the newton onset slope just above the first threshold + finite = np.sort(thresholds[thresholds < np.inf]) + d0 = thr0 + 0.4 * ((finite[1] - thr0) if finite.size > 1 else 0.3 * thr0) + df = compute_modal_intensities_full_salt_newton(g, tdf.copy(), d0, D0_steps=4) + cols = sorted( + c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities" + ) + a = np.nan_to_num(df.loc[t0, [("modal_intensities", c) for c in cols]].to_numpy(float)) + newton_slope = a[-1] / (cols[-1] - thr0) + assert 0.8 < newton_slope / linear_slope < 1.2 + + def test_two_mode_competition_negative_kink(self): + """Ge-Chong-Stone two-mode regression (PRA 82, 063824, Fig. 6): when the + second mode crosses its *interacting* threshold the operator solve must + (a) lase it -- the active set is found self-consistently, with the + within-edge hole burning resolved so the clamped background does not + freeze candidates below threshold -- and (b) suppress the dominant mode + below its un-kinked single-mode line (the negative competition kink). + Guards the active-set regression where a simultaneous multi-mode add + collapsed to single-mode lasing and the dominant rode its un-kinked line + (the full suite stayed green through that regression -- this test is the + pin).""" + from netsalt.modes import ( + compute_modal_intensities, + compute_modal_intensities_full_salt_newton, + compute_mode_competition_matrix, + ) + + g, tdf = _independent_lasing_fixture() + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + t0 = int(np.argmin(thresholds)) + thr0 = float(thresholds[t0]) + + # locate the second mode's interacting threshold from the linear/SPA sweep + T = compute_mode_competition_matrix(g, tdf) + d0_second = None + for d0 in np.linspace(thr0, 4.0 * thr0, 40): + lin = compute_modal_intensities(tdf.copy(), float(d0), T) + cols = [c for c in lin.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + last = np.nan_to_num(lin[max(cols, key=lambda c: c[1])].to_numpy(float)) + if np.sum(last > 1e-12) >= 2: + d0_second = float(d0) + break + assert d0_second is not None, "fixture must lase >=2 modes in the linear model" + + d0_max = 1.25 * d0_second + df = compute_modal_intensities_full_salt_newton(g, tdf.copy(), d0_max, D0_steps=8) + cols = sorted( + c[1] for c in df.columns if isinstance(c, tuple) and c[0] == "modal_intensities" + ) + curves = np.nan_to_num(df[[("modal_intensities", c) for c in cols]].to_numpy(float)) + # (a) the second mode lases + assert np.sum(curves[:, -1] > 1e-3 * curves[:, -1].max()) >= 2, ( + "operator solve must lase the second mode past its interacting threshold" + ) + # (b) negative kink: the dominant ends below its un-kinked single-mode + # extrapolation (onset slope from the first above-threshold points) + dom = curves[t0] + on = np.where(dom > 0)[0] + assert len(on) >= 3 + i0, i1 = on[0], on[1] + onset_slope = (dom[i1] - dom[i0]) / (cols[i1] - cols[i0]) + unkinked = dom[i0] + onset_slope * (cols[-1] - cols[i0]) + assert dom[-1] < 0.98 * unkinked, ( + "dominant mode must be suppressed below its un-kinked line " + f"(got {dom[-1]:.4g} vs un-kinked {unkinked:.4g})" + ) + + def test_linear_keeps_mode_ordering(self): + """The lowest-threshold (dominant) mode stays dominant far above + threshold in the linear model.""" + from netsalt.modes import compute_modal_intensities, compute_mode_competition_matrix + + g, tdf = _independent_lasing_fixture() + thresholds = np.asarray(tdf["lasing_thresholds"]).ravel() + assert np.sum(thresholds < np.inf) >= 2, "need >=2 lasing modes to test ordering" + t0 = int(np.argmin(thresholds)) + d0_max = 2.5 * float(thresholds[t0]) + + T = compute_mode_competition_matrix(g, tdf) + lin = compute_modal_intensities(tdf.copy(), d0_max, T) + cols = [c for c in lin.columns if isinstance(c, tuple) and c[0] == "modal_intensities"] + last = np.nan_to_num(lin[max(cols, key=lambda c: c[1])].to_numpy(float)) + assert int(np.argmax(last)) == t0