[eme] characterize orientation bias under modal truncation

Compare mixed and Galerkin testing against exact polarization/Fresnel interfaces and independent driven FDFD solves. Show complete-basis recovery, omitted-channel power, and the limits of longitudinal refinement, orientation averaging, and Euclidean least squares. Keep the default formulation unchanged and document its approximation limits.

Validation: 1689 tests passed, 584 skipped; Ruff and documentation build passed.

README.md

@@ -49,7 +49,7 @@ linear systems, ideally with double precision.
 
 **Requirements:**
 
-* python >=3.11
+* python >=3.12
 * numpy
 * scipy
 
@@ -118,8 +118,12 @@ python3 -m pytest -q | tee test_results.txt
 - `meanas.fdfd`: frequency-domain wave equations, sparse operators, SCPML, and
   iterative solves for driven problems.
 - `meanas.fdfd.waveguide_2d` / `meanas.fdfd.waveguide_3d`: waveguide mode
-  solvers, mode-source construction, and overlap windows for port-based
+  solvers, including targeted complex bases from the default half-sized
-  excitation and analysis.
+  beta-squared or explicit signed-beta formulations, transpose-left and
+  Lorentz-dual projection fields on a shared field-array axis, separately
+  returned PML-mode diagnostics, overlap-based tracking,
+  mode-source construction, and overlap windows for port-based excitation and
+  analysis.
 - `meanas.fdtd`: Yee-step updates, CPML boundaries, flux/energy accounting, and
   on-the-fly phasor extraction for comparing time-domain runs against FDFD.
 - `meanas.fdmath`: low-level finite-difference operators, vectorization helpers,
@@ -172,8 +176,10 @@ The tracked examples under `examples/` are the intended entry points for users:
 - `examples/eme_bend.py`: straight-to-bent waveguide mode matching with
   cylindrical bend modes, interface scattering, and a cascaded bend-network
   example built with `scikit-rf`.
-- `examples/eme_taper_cmt.py`: sampled cross-section local-mode CMT for a
+- `examples/eme_taper_cmt.py`: sampled-interface EME with stable S cascading for a
   continuously varying rib-waveguide taper.
+- `examples/eme_taper_pml_modes.py`: targeted complex/PML mode solving,
+  tracking, and construction of EME taper sections.
 - `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap /
   Poynting analysis without a time-domain run.
 

docs/api/fdfd.md

@@ -13,3 +13,24 @@
 ::: meanas.fdfd.scpml
 
 ::: meanas.fdfd.farfield
+
+## SCPML conventions
+
+`stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the
+high-index face, matching CPML and propagation conventions. It returns copied
+widths by default; pass `copy=False` to request mutation explicitly.
+
+Both SCPML builders accept `kappa_function(depth)` for real stretching and
+`max_imag_ratio` to limit the final dimensionless `abs(Im(s))/Re(s)`. Start a
+mode-solver convergence study with `max_imag_ratio=1.0`: overly strong imaginary
+stretching can produce above-light-line PML artifacts. The cap is not a general
+spectral guarantee, and weaker absorption can increase boundary reflections.
+Check domain size, PML thickness, strength, mode localization, and residuals.
+The general scattering profile remains uncapped by default.
+
+`uniform_grid_scpml` builds unit widths by applying `stretch_with_scpml` to
+both faces of each enabled axis. Both builders divide imaginary grading by
+physical layer thickness and sample electric widths at cell centers and
+magnetic widths at cell edges. This corrects the uniform builder's previous
+thickness-independent strength and reversed E/H sample offsets; retune existing
+uniform-builder profiles using domain/PML convergence checks.

docs/api/fdtd.md

@@ -13,3 +13,12 @@
 ::: meanas.fdtd.energy
 
 ::: meanas.fdtd.phasor
+
+CPML derivatives use electric widths for forward differences and magnetic
+widths for backward differences. `cpml_params(..., dx=1.0)` specifies the uniform
+physical cell width normal to a face; arbitrary nonuniform PML grading is not
+implemented by this constructor. Scale the timestep consistently with the grid.
+
+`ln_R_per_layer` must be finite and strictly negative. Zero attenuation raises
+`ValueError`; disable individual CPML faces with `None` in the parameter table
+passed to `updates_with_cpml`.

docs/api/waveguides.md

@@ -5,3 +5,149 @@
 ::: meanas.fdfd.waveguide_3d
 
 ::: meanas.fdfd.waveguide_cyl
+
+## Overlap measures and normalization
+
+`waveguide_2d.inner_product(..., measure="physical", conj_h=True)` uses real
+Yee widths for physical Poynting power. `measure="stretched"` keeps complex
+widths for coordinate reaction. `lorentz_metric` and `lorentz_overlap` default
+to the stretched measure. Use the same grid and measure when constructing and
+using dual fields. A real overlap grid still defines a consistent constructed
+dual basis, even when the eigenproblem uses a separate PML grid.
+
+`normalized_fields_e/h(..., normalization="lorentz")` normalizes unconjugated
+self-reaction to one. This fixes phase up to a sign and is distinct from unit
+forward power, which remains the default. Zero-power or self-orthogonal modes
+raise `ValueError`; open self-orthogonal modes require paired primal/dual
+normalization through `ModeBasis`.
+
+## Targeted and reproducible mode solving
+
+For modes near an effective index guess, use `solve_modes_near` with
+`target_wavenumbers=[omega * n_guess]`. It targets the electric beta-squared
+operator at `(omega * n_guess)**2`; `solve_modes_near_beta` instead targets the
+first-order signed-beta operator. Both use the LU shift-invert backend and
+support `eigensolver_opts={"v0": starting_vector}`. Factorization cost and
+memory depend on the grid, so this is not universally faster than indexed
+solving. `solve_modes` keeps its global mode-number semantics and now accepts
+`rng` or `v0` for reproducible initialization.
+
+`waveguide_3d.mode_from_2d` accepts a supplied transverse electric eigenvector
+and reuses field reconstruction, normalization, numerical-dispersion
+correction, Yee staggering, rotation, and full-grid embedding without solving
+an eigenproblem.
+
+## Interface and sampled-taper scattering
+
+All scattering matrices use outgoing rows and incoming columns, ordered
+`[left, right]`:
+
+```text
+[b_left ]   [R_left        T_right_to_left] [a_left ]
+[b_right] = [T_left_to_right       R_right] [a_right]
+```
+
+`get_s` and `get_abcd` require explicit `dxes`. Mixed Petrov-Galerkin testing is
+the default: left electric duals test H continuity, and right magnetic duals
+test E continuity. Both incidence directions are solved together with full
+projection matrices. Mixed scattering supports unequal mode counts.
+`formulation="galerkin"` selects equal-count left testing of both equations.
+Neither formulation silently clips transmission. The optional `force_nogain`
+and `force_reciprocal` flags are numerical projections, not general physical
+tests for arbitrary open or evanescent modal coordinates.
+
+`get_taper_s` replaces finite-difference CMT with half propagation in the left
+basis, a midpoint interface, and half propagation in the right basis for each
+sampled interval. Redheffer composition avoids growing passive backward
+transfer branches. Refine section sampling and modal/PML bases to establish
+convergence. The ineffective `rtol`, `atol`, and `max_step` arguments have been
+removed. The existing `eme_taper_cmt.py` example filename is retained.
+
+```python
+prepared = eme.prepare_taper(sections, dxes)
+s_original = eme.evaluate_taper_s(prepared)
+s_longer = eme.evaluate_taper_s(prepared, z_coords=prepared.z_coords * 1.2)
+s_changed_beta = eme.evaluate_taper_s(
+    prepared, wavenumber_matrices=beta_matrices,
+)
+```
+
+`PreparedTaper` owns read-only copies of phase-aligned interface data. Beta
+matrix overrides have shape `(section, mode, mode)` in the original supplied
+section bases; recorded basis transforms are applied automatically. Explicit
+backward branches require backward beta overrides when forward beta is changed.
+Evaluation performs no field overlaps or eigensolves. Length sweeps reuse the
+sampled model; wavelength sweeps with reused fields are explicitly frozen-basis
+approximations. Reprepare when fields, duals, grids, or formulations change.
+
+`cascade_s(left, right, internal_modes=...)` connects spatially ordered S
+matrices, including unequal external ports. Public `abcd_to_s` and `s_to_abcd`
+use `[f_right, b_right] = ABCD @ [f_left, b_left]` and require equal port counts.
+`get_taper_abcd` converts the stable S result; zero transmission or an
+unrepresentable transfer matrix raises `LinAlgError`. Prefer S for long lossy
+devices. Singular interface and cascade systems also raise contextual errors
+rather than silently substituting pseudoinverses.
+
+::: meanas.fdfd.eme
+
+## Orientation bias under modal truncation
+
+Mixed testing is a projected approximation: H continuity is tested on the left
+and E continuity on the right. Swapping the ports changes the testing spaces
+when the retained modal subspaces are incomplete. This effect survives checks
+of reciprocity and unitarity.
+
+A reproducible Maxwell counterexample is a uniform isotropic medium with two
+degenerate linear polarizations. Keep one polarization on each port, rotated
+so their overlap is `a=0.8`. There is no material interface: true reflection is
+zero and transmission into the retained receiver is `a`. The other polarization
+carries the remaining transmitted power, `1-a**2=0.36`.
+
+| Calculation | R_left | R_right | T_left_to_right | T_right_to_left |
+| --- | ---: | ---: | ---: | ---: |
+| Exact full basis, then restrict ports | 0 | 0 | 0.8 | 0.8 |
+| Truncated mixed | 0.219512 | -0.219512 | 0.975610 | 0.975610 |
+| Swap ports, mixed, then reorder | -0.219512 | 0.219512 | 0.975610 | 0.975610 |
+| Truncated left Galerkin | 0 | 0 | 1.25 | 0.8 |
+| Euclidean E/H least squares | 0 | 0 | 0.8 | 0.8 |
+
+An independent driven 3D FDFD solve with PML agrees with the exact reference
+from either incidence direction: transmission is 0.8 within `1e-12`, with
+reflection below `2e-8` in the test fixture. Keeping both polarizations makes
+mixed testing match the exact complete scattering matrix and restores
+port-relabeling covariance.
+
+For the single retained polarization, mixed testing yields
+`r=(1-a**2)/(1+a**2)` and `t=2*a/(1+a**2)`. The spurious reflection is quadratic
+in a small angle between adjacent retained polarizations. However, longitudinal
+refinement alone is insufficient: rotating a single retained vector through a
+uniform degenerate medium makes the sampled chain approach artificial following
+of that vector (`|T|` approaches 1), while the physical endpoint overlap remains
+0.8. Keeping the complete degenerate block recovers the physical result at every
+sampling density tested.
+
+Two tempting corrections are insufficient. Averaging the two port orientations
+cancels reflection in this example but leaves transmission at 0.975610. Plain
+Euclidean E/H least squares gives the correct uniform-medium result, but is not
+a general replacement. A second fixture uses exact degenerate modes of diagonal
+positive epsilon/mu tensors, with electric scaling `(1,2)` on the left and
+`(1,3)` on the right and a 30-degree polarization rotation on both ports. After
+single-mode truncation, least squares has unequal directional transmissions
+(approximately 0.850839 and 1.092163) and a maximum singular value of 1.104947.
+The complete Fresnel solution is reciprocal and lossless.
+
+The mixed default remains unchanged. Retain complete coupled/degenerate blocks
+when the geometry can mix them, check modal-basis convergence as well as section
+sampling, and compare against independent full-field results. Neither orientation
+averaging nor gain/reciprocity postprocessing recovers omitted-channel physics.
+Orientation discrepancy is a useful diagnostic, not a bound on transmission
+error. These examples use real power-orthonormal propagating modes; for general
+complex modes, the reciprocity constraint involves the modal orthogonality
+matrix, as derived by [Svendsen, Haakestad, and Skaar](https://arxiv.org/abs/1301.2458).
+
+Reproduce the analytic, full-basis, candidate-formulation, taper-refinement, and
+driven-FDFD checks with:
+
+```sh
+OPENBLAS_NUM_THREADS=1 python -m pytest meanas/test/test_eme_orientation.py -q
+```

docs/index.md

@@ -28,7 +28,7 @@ Relevant starting examples:
   scattering between two nearby waveguide cross-sections
 - `examples/eme_bend.py` for straight-to-bent mode matching with cylindrical
   bend modes and a cascaded bend-network example
-- `examples/eme_taper_cmt.py` for local-mode CMT through sampled continuously
+- `examples/eme_taper_cmt.py` for sampled-interface EME through continuously
   varying taper cross-sections
 - `examples/fdfd.py` for direct frequency-domain waveguide excitation
 

examples/eme_bend.py

@@ -264,7 +264,7 @@ def main() -> None:
     skrf = require_optional('skrf', package_name='scikit-rf')
 
     grid, epsilon, dxes_2d, omega = build_geometry()
-    epsilon_slice = epsilon[:, :, :, 2]
+    epsilon_slice = epsilon[:, :, :, 1]
     rmin = RADIUS + grid.xyz[0].min()
 
     straight_modes, straight_wavenumbers = solve_straight_modes(epsilon_slice, omega=omega, dxes_2d=dxes_2d)

examples/eme_taper_cmt.py

@@ -1,5 +1,5 @@
 """
-Local-mode CMT example for a continuously varying rib-waveguide taper.
+Sampled-interface EME example for a continuously varying rib-waveguide taper.
 
 This example keeps geometry construction outside `meanas.fdfd.eme`: it samples a
 width taper at several cross-sections, solves and normalizes each local mode with
@@ -112,14 +112,15 @@ def print_summary(
     print('sampled taper effective indices:', ', '.join(f'{value:.5f}' for value in neffs))
     print(f'abrupt endpoint reflection       |S_00|^2          = {abs(abrupt_s[0, 0]) ** 2:.6f}')
     print(f'abrupt endpoint transmission     |S_{n_modes},0|^2 = {abs(abrupt_s[n_modes, 0]) ** 2:.6f}')
-    print(f'taper CMT reflection             |S_00|^2          = {abs(taper_s[0, 0]) ** 2:.6f}')
+    print(f'taper reflection             |S_00|^2          = {abs(taper_s[0, 0]) ** 2:.6f}')
-    print(f'taper CMT transmission           |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}')
+    print(f'taper transmission           |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}')
-    print(f'taper CMT total output power                       = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}')
+    print(f'taper total output power                       = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}')
 
 
 def main() -> None:
     sections, neffs, _omega, dxes_2d = solve_taper_sections()
-    taper_s = eme.get_taper_s(sections, dxes=dxes_2d)
+    prepared = eme.prepare_taper(sections, dxes=dxes_2d)
+    taper_s = eme.evaluate_taper_s(prepared)
     abrupt_s = eme.get_s(
         sections[0].modes,
         sections[0].wavenumbers,

examples/eme_taper_pml_modes.py

@@ -0,0 +1,110 @@
+"""Build tracked complex/PML modal bases for a slowly varying taper."""
+
+from __future__ import annotations
+
+import numpy
+
+from meanas.fdmath import vec
+from meanas.fdfd import eme, scpml, waveguide_2d
+
+
+WAVELENGTH = 1.55
+OMEGA = 2 * numpy.pi / WAVELENGTH
+SHAPE = (17, 17)
+SPACING = 0.25
+PML_CELLS = 3
+EPS_CLADDING = 1.45 ** 2
+EPS_CORE = 2.0 ** 2
+
+
+def build_grids() -> tuple[list[list[numpy.ndarray]], list[list[numpy.ndarray]]]:
+    """Build matching physical and complex-stretched transverse Yee grids.
+
+    Returns:
+        The real overlap grid followed by the PML-stretched solve grid.
+    """
+    physical = [[numpy.full(size, SPACING) for size in SHAPE] for _ in range(2)]
+    stretched = [[values.copy() for values in family] for family in physical]
+    for axis in (0, 1):
+        for polarity in (-1, +1):
+            stretched = scpml.stretch_with_scpml(
+                stretched,
+                axis=axis,
+                polarity=polarity,
+                omega=OMEGA,
+                epsilon_effective=EPS_CLADDING,
+                thickness=PML_CELLS,
+                max_imag_ratio=1.0,
+            )
+    return physical, stretched
+
+
+def build_epsilon(edge_fill: float) -> numpy.ndarray:
+    """Construct and vectorize one taper cross-section permittivity profile.
+
+    Args:
+        edge_fill: Fractional core fill used for the two cells beside the
+            central three-cell-wide core.
+
+    Returns:
+        Vectorized three-component relative permittivity field.
+    """
+    epsilon = numpy.full((3, *SHAPE), EPS_CLADDING)
+    center_x, center_y = (size // 2 for size in SHAPE)
+    epsilon[:, center_x - 1:center_x + 2, center_y] = EPS_CORE
+    mixed_epsilon = EPS_CLADDING + edge_fill * (EPS_CORE - EPS_CLADDING)
+    epsilon[:, center_x - 1:center_x + 2, center_y - 1] = mixed_epsilon
+    epsilon[:, center_x - 1:center_x + 2, center_y + 1] = mixed_epsilon
+    return vec(epsilon)
+
+
+def main() -> None:
+    """Solve, track, and summarize the sampled complex modal bases."""
+    physical_dxes, solve_dxes = build_grids()
+    z_coords = numpy.array([0.0, 5.0, 10.0])
+    epsilon_samples = [build_epsilon(edge_fill) for edge_fill in (0.30, 0.35, 0.40)]
+    solved = [
+        waveguide_2d.solve_modes_near(
+            OMEGA,
+            solve_dxes,
+            epsilon,
+            target_wavenumbers=[1.45 * OMEGA],
+            modes_per_target=3,
+            candidate_margin=3,
+            overlap_dxes=physical_dxes,
+        )
+        for epsilon in epsilon_samples
+    ]
+    bases = [basis for basis, _diagnostics in solved]
+    tracked = waveguide_2d.track_mode_path(bases, overlap_dxes=physical_dxes)
+    tracked_diagnostics = [
+        waveguide_2d.compute_mode_diagnostics(
+            basis,
+            OMEGA,
+            solve_dxes,
+            epsilon,
+            target_wavenumbers=[1.45 * OMEGA],
+            overlap_dxes=physical_dxes,
+        )
+        for basis, epsilon in zip(tracked, epsilon_samples, strict=True)
+    ]
+    sections = eme.mode_sections_from_bases(z_coords, tracked)
+    prepared = eme.prepare_taper(sections, physical_dxes)
+    # Fields and interface overlaps are reused: these are frozen-basis length
+    # sweeps. Wavelength-dependent beta overrides would also freeze the fields.
+    length_sweep = [eme.evaluate_taper_s(prepared, z_coords=z_coords * scale)
+                    for scale in (0.8, 1.0, 1.2)]
+
+    for z_coord, basis, diagnostics in zip(z_coords, tracked, tracked_diagnostics, strict=True):
+        print(
+            f'z={z_coord:5.1f}: '
+            f'beta={basis.wavenumbers[0]:.6g}, '
+            f'PML participation={diagnostics.pml_participation[0]:.3f}, '
+            f'right residual={diagnostics.right_residuals[0]:.2e}',
+        )
+    print(f'constructed {len(sections)} tracked EME sections with {len(sections[0].modes)} modes each')
+    print('frozen-basis length sweep |T_00|:', [float(abs(ss[len(tracked[0]), 0])) for ss in length_sweep])
+
+
+if __name__ == '__main__':
+    main()

examples/fdtd.py

@@ -145,7 +145,7 @@ def main() -> None:
 
     # PMLs in every direction
     pml_params = [
-        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_air ** 2)
+        [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, thickness=pml_thickness, epsilon_eff=n_air ** 2)
          for pp in (-1, +1)]
         for dd in range(3)]
     update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon, dtype=complex)
@@ -203,7 +203,7 @@ def main() -> None:
     dxes_fdfd = copy.deepcopy(dxes)
     for pp in (-1, +1):
         for dd in range(3):
-            stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_air ** 2, thickness=pml_thickness)
+            dxes_fdfd = stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_air ** 2, thickness=pml_thickness)
     # Compare the extracted phasor to the FDFD operator on the stretched grid,
     # not the unstretched Yee spacings used by the raw time-domain update.
     A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon)

examples/waveguide.py

@@ -273,7 +273,7 @@ def main2():
 
     # PMLs in every direction
     pml_params = [
-        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_cladding ** 2)
+        [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, thickness=pml_thickness, epsilon_eff=n_cladding ** 2)
          for pp in (-1, +1)]
         for dd in range(3)]
     update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon, dtype=complex)
@@ -330,7 +330,7 @@ def main2():
     dxes_fdfd = copy.deepcopy(dxes)
     for pp in (-1, +1):
         for dd in range(3):
-            stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_cladding ** 2, thickness=pml_thickness)
+            dxes_fdfd = stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_cladding ** 2, thickness=pml_thickness)
     # Residuals must be checked on the stretched FDFD grid, because the FDTD run
     # already includes those same absorbing layers through CPML.
     A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon)

meanas/eigensolvers.py

@@ -2,17 +2,153 @@
 Solvers for eigenvalue / eigenvector problems
 """
 from collections.abc import Callable
+import dataclasses
 import numpy
 from numpy.typing import NDArray, ArrayLike
 from numpy.linalg import norm
 from scipy import sparse
 import scipy.sparse.linalg as spalg
+from scipy.optimize import linear_sum_assignment
+
+
+@dataclasses.dataclass(frozen=True, eq=False)
+class TargetedEigensolution:
+    """Paired right and transpose-left eigenvectors near a spectral target.
+
+    Attributes:
+        eigenvalues: Eigenvalues associated with the columns of both vector
+            matrices.
+        right_eigenvectors: Right eigenvectors stored column-wise.
+        left_eigenvectors: Bilinear-left eigenvectors stored column-wise. These
+            satisfy a transpose, rather than Hermitian-transpose, eigenproblem.
+    """
+
+    eigenvalues: NDArray[numpy.complex128]
+    right_eigenvectors: NDArray[numpy.complex128] = dataclasses.field(repr=False)
+    left_eigenvectors: NDArray[numpy.complex128] = dataclasses.field(repr=False)
+
+
+type TargetedEigensolver = Callable[..., TargetedEigensolution]
+
+
+def scipy_eigs_near(
+        operator: sparse.sparray,
+        target: complex,
+        how_many: int,
+        *,
+        tolerance: float = 1e-10,
+        max_iterations: int | None = None,
+        ncv: int | None = None,
+        pairing_rtol: float = 1e-7,
+        v0: ArrayLike | None = None,
+        ) -> TargetedEigensolution:
+    """
+    Solve for paired right and bilinear-left eigenvectors near ``target``.
+
+    One sparse LU factorization of ``operator - target * I`` is reused for the
+    right solve and for the transpose-left solve.  The left vectors therefore
+    satisfy ``left.T @ operator = eigenvalue * left.T``; this deliberately uses
+    a transpose rather than a Hermitian transpose so the result is suitable for
+    unconjugated Lorentz-reciprocity projections.
+
+    Args:
+        operator: Square sparse operator to diagonalize.
+        target: Shift-invert spectral target.
+        how_many: Number of eigenpairs to return.
+        tolerance: ARPACK convergence tolerance for both solves.
+        max_iterations: Optional ARPACK iteration limit.
+        ncv: Optional number of Arnoldi vectors. A size-aware default is used
+            when omitted.
+        pairing_rtol: Relative tolerance for matching right and transpose-left
+            eigenvalues.
+        v0: Optional finite starting vector shared by both solves.
+
+    Returns:
+        Paired eigenvalues and column-wise right and transpose-left vectors.
+
+    Raises:
+        ValueError: If the dimensions, counts, tolerances, or starting vector
+            are invalid.
+        RuntimeError: If the independently solved spectra cannot be paired.
+    """
+    if how_many < 1:
+        raise ValueError('how_many must be positive')
+    if not numpy.isfinite(tolerance) or tolerance <= 0:
+        raise ValueError('tolerance must be finite and positive')
+    if not numpy.isfinite(pairing_rtol) or pairing_rtol <= 0:
+        raise ValueError('pairing_rtol must be finite and positive')
+
+    operator = sparse.csc_array(operator, dtype=complex)
+    size = operator.shape[0]
+    if operator.shape != (size, size):
+        raise ValueError('operator must be square')
+    if how_many >= size - 1:
+        raise ValueError('how_many must be smaller than operator dimension minus one')
+
+    shifted = operator - target * sparse.eye_array(size, format='csc', dtype=complex)
+    factor = spalg.splu(shifted)
+    dtype = numpy.dtype(numpy.complex128)
+    inverse = spalg.LinearOperator(
+        operator.shape,
+        matvec=lambda vector: factor.solve(vector, trans='N'),
+        dtype=dtype,
+    )
+    inverse_transpose = spalg.LinearOperator(
+        operator.shape,
+        matvec=lambda vector: factor.solve(vector, trans='T'),
+        dtype=dtype,
+    )
+    if ncv is None:
+        ncv = min(size, max(2 * how_many + 1, 4 * how_many))
+    start = None if v0 is None else numpy.asarray(v0, dtype=complex)
+    if start is not None and (start.shape != (size,) or not numpy.isfinite(start).all()):
+        raise ValueError('v0 must be a finite vector matching the operator dimension')
+
+    eigenvalues, right = spalg.eigs(
+        operator,
+        k=how_many,
+        sigma=target,
+        which='LM',
+        OPinv=inverse,
+        tol=tolerance,
+        maxiter=max_iterations,
+        ncv=ncv,
+        v0=start,
+    )
+    left_values, left = spalg.eigs(
+        operator.T,
+        k=how_many,
+        sigma=target,
+        which='LM',
+        OPinv=inverse_transpose,
+        tol=tolerance,
+        maxiter=max_iterations,
+        ncv=ncv,
+        v0=start,
+    )
+
+    rows, columns = linear_sum_assignment(numpy.abs(eigenvalues[:, None] - left_values[None, :]))
+    paired_columns = columns[numpy.argsort(rows)]
+    left_values = left_values[paired_columns]
+    left = left[:, paired_columns]
+    pairing_error = numpy.abs(eigenvalues - left_values)
+    pairing_scale = numpy.maximum(1.0, numpy.maximum(numpy.abs(eigenvalues), numpy.abs(left_values)))
+    if numpy.any(pairing_error > pairing_rtol * pairing_scale):
+        raise RuntimeError('could not pair right and transpose-left eigenvalues')
+
+    return TargetedEigensolution(
+        eigenvalues=numpy.asarray(eigenvalues, dtype=complex),
+        right_eigenvectors=numpy.asarray(right, dtype=complex),
+        left_eigenvectors=numpy.asarray(left, dtype=complex),
+    )
 
 
 def power_iteration(
         operator: sparse.spmatrix,
         guess_vector: NDArray[numpy.complex128] | None = None,
         iterations: int = 20,
+        *,
+        rng: numpy.random.Generator | int | None = None,
         ) -> tuple[complex, NDArray[numpy.complex128]]:
     """
     Use power iteration to estimate the dominant eigenvector of a matrix.
@@ -21,15 +157,18 @@ def power_iteration(
        operator: Matrix to analyze.
        guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
        iterations: Number of iterations to perform. Default 20.
+       rng: Generator or seed used only when ``guess_vector`` is omitted.
 
     Returns:
         (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
     """
-    rng = numpy.random.default_rng()
     if guess_vector is None:
+        rng = numpy.random.default_rng(rng)
         v = rng.random(operator.shape[0]) + 1j * rng.random(operator.shape[0])
     else:
-        v = guess_vector
+        v = numpy.array(guess_vector, dtype=complex, copy=True)
+    if v.shape != (operator.shape[0],) or not numpy.isfinite(v).all() or not numpy.any(v):
+        raise ValueError('guess_vector must be finite, nonzero, and match the operator dimension')
 
     for _ in range(iterations):
         v = operator @ v
@@ -67,12 +206,14 @@ def rayleigh_quotient_iteration(
         (operator - sparse.eye_array(operator.shape[0]))
 
         def shift(eigval: float) -> sparse.sparray:
+            """Construct a sparse scalar identity shift."""
             return eigval * sparse.eye_array(operator.shape[0])
 
         if solver is None:
             solver = spalg.spsolve
     except TypeError:
         def shift(eigval: float) -> spalg.LinearOperator:
+            """Construct a matrix-free scalar identity shift."""
             return spalg.LinearOperator(
                 shape=operator.shape,
                 dtype=operator.dtype,
@@ -80,6 +221,7 @@ def rayleigh_quotient_iteration(
                 )
         if solver is None:
             def solver(A: spalg.LinearOperator, b: ArrayLike) -> NDArray[numpy.complex128]:
+                """Solve a matrix-free shifted system with BiCGSTAB."""
                 return spalg.bicgstab(A, b)[0]
     assert solver is not None
 
@@ -100,6 +242,9 @@ def signed_eigensolve(
         operator: sparse.spmatrix | spalg.LinearOperator,
         how_many: int,
         negative: bool = False,
+        *,
+        rng: numpy.random.Generator | int | None = None,
+        v0: ArrayLike | None = None,
         ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
     """
     Find the largest-magnitude positive-only (or negative-only) eigenvalues and
@@ -110,13 +255,24 @@ def signed_eigensolve(
         how_many: How many eigenvalues to find.
         negative: Whether to find negative-only eigenvalues.
                   Default False (positive only).
+        rng: Optional generator or seed for reproducible initialization.
+        v0: Finite nonzero starting vector, taking precedence over ``rng``.
 
     Returns:
         (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
         `eigenvectors[:, k]` corresponds to the k-th eigenvalue
     """
     # Use power iteration to estimate the dominant eigenvector
-    lm_eigval, _ = power_iteration(operator)
+    if v0 is None:
+        generator = numpy.random.default_rng(rng)
+        start = generator.random(operator.shape[0]) + 1j * generator.random(operator.shape[0])
+    else:
+        start = numpy.array(v0, dtype=complex, copy=True)
+    lm_eigval, _ = power_iteration(operator, guess_vector=start)
+    if not numpy.issubdtype(operator.dtype, numpy.complexfloating):
+        start = start.real
+        if not numpy.any(start):
+            raise ValueError('v0 must have a nonzero real part for a real operator')
 
     '''
     Shift by the absolute value of the largest eigenvalue, then find a few of the
@@ -135,9 +291,11 @@ def signed_eigensolve(
         shifted_operator = operator + spalg.LinearOperator(shape=operator.shape,
                                                            matvec=lambda v: shift * v)
 
-    shifted_eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=2 * how_many + 50)
+    shifted_eigenvalues, eigenvectors = spalg.eigs(
+        shifted_operator, which='LM', k=how_many,
+        ncv=min(operator.shape[0], 2 * how_many + 50), v0=start,
+    )
     eigenvalues = shifted_eigenvalues - shift
 
     k = eigenvalues.argsort()
     return eigenvalues[k], eigenvectors[:, k]
-

meanas/fdfd/scpml.py

@@ -38,132 +38,148 @@ def prepare_s_function(
     return s_factor
 
 
+def _stretch_factors(
+        distance: NDArray,
+        correction: float,
+        s_function: s_function_t,
+        kappa_function: s_function_t | None,
+        max_imag_ratio: float | None,
+        ) -> NDArray[numpy.complex128]:
+    """Evaluate the final dimensionless coordinate stretch."""
+    kappa = numpy.asarray(1.0 if kappa_function is None else kappa_function(distance))
+    sigma = numpy.asarray(s_function(distance)) / correction
+    if numpy.iscomplexobj(kappa) or not numpy.isfinite(kappa).all() or (kappa < 1).any():
+        raise ValueError('kappa_function must return finite real stretches at least one')
+    if numpy.iscomplexobj(sigma) or not numpy.isfinite(sigma).all():
+        raise ValueError('s_function must return finite real values')
+    if max_imag_ratio is not None:
+        sigma = numpy.clip(sigma, -max_imag_ratio * kappa, max_imag_ratio * kappa)
+    return kappa + 1j * sigma
+
+
+def _validate_profile(omega: float, epsilon_effective: float, max_imag_ratio: float | None) -> None:
+    if not numpy.isfinite(omega) or not numpy.isreal(omega) or omega <= 0:
+        raise ValueError('omega must be finite, real, and positive')
+    if not numpy.isfinite(epsilon_effective) or epsilon_effective <= 0:
+        raise ValueError('epsilon_effective must be finite and positive')
+    if max_imag_ratio is not None and (not numpy.isfinite(max_imag_ratio) or max_imag_ratio <= 0):
+        raise ValueError('max_imag_ratio must be finite and positive')
+
+
 def uniform_grid_scpml(
         shape: Sequence[int],
         thicknesses: Sequence[int],
         omega: float,
         epsilon_effective: float = 1.0,
         s_function: s_function_t | None = None,
-        ) -> list[list[NDArray[numpy.float64]]]:
+        *,
-    """
+        kappa_function: s_function_t | None = None,
-    Create dx arrays for a uniform grid with a cell width of 1 and a pml.
+        max_imag_ratio: float | None = None,
+        ) -> list[list[NDArray]]:
+    """Create unit-cell Yee widths with SCPML on both faces of each axis.
 
-    If you want something more fine-grained, check out `stretch_with_scpml(...)`.
+    ``shape`` and ``thicknesses`` contain three axis sizes and per-face PML
+    thicknesses; zero disables an axis. ``s_function`` specifies the imaginary
+    grading, defaulting to prepare_s_function(). ``kappa_function`` optionally
+    returns dimensionless real stretch values >= 1 at normalized depths.
+    Each face is built by stretch_with_scpml, including its physical-thickness
+    scaling and electric-center / magnetic-edge sampling convention.
 
-    Args:
+    ``max_imag_ratio`` optionally bounds abs(Im(s))/Re(s) after frequency
-        shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick)
+    scaling. A value of 1 is a useful mode-solver starting point, not a proof
-        thicknesses: `[th_x, th_y, th_z]`
+    against spectral artifacts. Converge domain size, thickness, and strength;
-                     Thickness of the PML in each direction.
+    reducing absorption can increase boundary reflections. The general
-                     Both polarities are added.
+    scattering profile remains uncapped by default.
-                     Each th_ of pml is applied twice, once on each edge of the grid along the given axis.
-                     `th_*` may be zero, in which case no pml is added.
-        omega: Angular frequency for the simulation
-        epsilon_effective: Effective epsilon of the PML. Match this to the material
-                            at the edge of your grid.
-                            Default 1.
-        s_function: created by `prepare_s_function(...)`, allowing customization of pml parameters.
-                    Default uses `prepare_s_function()` with no parameters.
-
-    Returns:
-        Complex cell widths (dx_lists_mut) as discussed in `meanas.fdmath.types`.
     """
+    _validate_profile(omega, epsilon_effective, max_imag_ratio)
+    shape, thicknesses = tuple(shape), tuple(thicknesses)
+    if len(shape) != 3 or len(thicknesses) != 3:
+        raise ValueError('shape and thicknesses must contain three entries')
+    if any(not isinstance(n, int | numpy.integer) or isinstance(n, bool) or n <= 0 for n in shape):
+        raise ValueError('shape entries must be positive integers')
+    if any(not isinstance(t, int | numpy.integer) or isinstance(t, bool) or t < 0 or 2 * t > n
+           for n, t in zip(shape, thicknesses, strict=True)):
+        raise ValueError('PML thicknesses must be nonnegative integers fitting both faces')
     if s_function is None:
         s_function = prepare_s_function()
-
+    result = [[numpy.ones(n) for n in shape] for _ in range(2)]
-    shape = tuple(shape)
+    for axis, thickness in enumerate(thicknesses):
-    thicknesses = tuple(thicknesses)
+        if not thickness:
-
+            continue
-    # Normalized distance to nearest boundary
+        for polarity in (-1, 1):
-    def ll(u: NDArray[numpy.float64], n: int, t: int) -> NDArray[numpy.float64]:
+            stretch_with_scpml(
-        return ((t - u).clip(0) + (u - (n - t)).clip(0)) / t
+                result, axis, polarity, omega, epsilon_effective, thickness,
-
+                s_function, kappa_function=kappa_function,
-    dx_a = [numpy.array(numpy.inf)] * 3
+                max_imag_ratio=max_imag_ratio, copy=False,
-    dx_b = [numpy.array(numpy.inf)] * 3
+            )
-
+    return result
-    # divide by this to adjust for epsilon_effective and omega
-    s_correction = numpy.sqrt(epsilon_effective) * numpy.real(omega)
-
-    for k, th in enumerate(thicknesses):
-        s = shape[k]
-        if th > 0:
-            sr = numpy.arange(s, dtype=numpy.float64)
-            dx_a[k] = 1 + 1j * s_function(ll(sr,       s, th)) / s_correction
-            dx_b[k] = 1 + 1j * s_function(ll(sr + 0.5, s, th)) / s_correction
-        else:
-            dx_a[k] = numpy.ones((s,))
-            dx_b[k] = numpy.ones((s,))
-    return [dx_a, dx_b]
 
 
 def stretch_with_scpml(
-        dxes: list[list[NDArray[numpy.float64]]],
+        dxes: list[list[NDArray]],
         axis: int,
         polarity: int,
         omega: float,
         epsilon_effective: float = 1.0,
         thickness: int = 10,
         s_function: s_function_t | None = None,
-        ) -> list[list[NDArray[numpy.float64]]]:
+        *,
-    """
+        kappa_function: s_function_t | None = None,
-        Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis.
+        max_imag_ratio: float | None = None,
+        copy: bool = True,
+        ) -> list[list[NDArray]]:
+    """Add SCPML to one face, returning independent widths by default.
 
-        Args:
+    Args:
-            dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types`
+        dxes: Yee widths [dx_e, dx_h]; two or three spatial axes are supported.
-            axis: axis to stretch (0=x, 1=y, 2=z)
+        axis: Axis normal to the absorbing face.
-            polarity: direction to stretch (-1 for -ve, +1 for +ve)
+        polarity: -1 selects the low-index face; +1 selects the high-index face.
-            omega: Angular frequency for the simulation
+        omega: Positive real angular frequency.
-            epsilon_effective: Effective epsilon of the PML. Match this to the material at the
+        epsilon_effective: Positive effective permittivity at the boundary.
-                               edge of your grid. Default 1.
+        thickness: Number of cells on this face; zero is a no-op.
-            thickness: number of cells to use for pml (default 10)
+        s_function: Imaginary grading at normalized depths, from prepare_s_function.
-            s_function: Created by `prepare_s_function(...)`, allowing customization
+        kappa_function: Optional dimensionless real stretch profile, >= 1.
-                        of pml parameters. Default uses `prepare_s_function()` with no parameters.
+        max_imag_ratio: Optional final abs(Im(s))/Re(s) bound. A value of 1 is
+            a mode-solver starting configuration; check PML/domain convergence.
+        copy: Return copied widths by default. False explicitly mutates dxes.
 
-        Returns:
+    Coordinates and physical PML thickness are obtained from the real input
-            Complex cell widths (dx_lists_mut) as discussed in `meanas.fdmath.types`.
+    electric widths. Multiple nonoverlapping faces can be added successively.
-            Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed.
+    The imaginary sign follows the exp(+i omega t) convention used by meanas.
     """
+    _validate_profile(omega, epsilon_effective, max_imag_ratio)
+    if len(dxes) != 2 or len(dxes[0]) != len(dxes[1]) or axis not in range(len(dxes[0])):
+        raise ValueError('axis must select a spatial axis in both Yee width families')
+    if polarity not in (-1, 1):
+        raise ValueError('polarity must be -1 or +1')
+    if not isinstance(thickness, int | numpy.integer) or isinstance(thickness, bool):
+        raise TypeError('thickness must be an integer')
+    n = len(dxes[0][axis])
+    if thickness < 0 or thickness > n or len(dxes[1][axis]) != n:
+        raise ValueError('thickness must fit the matching Yee width arrays')
+    widths = [[values.copy() for values in family] for family in dxes] if copy else dxes
+    if thickness == 0:
+        return widths
     if s_function is None:
         s_function = prepare_s_function()
-
+    da, db = (numpy.array(family[axis], dtype=complex, copy=True) for family in widths)
-    dx_ai = dxes[0][axis].astype(complex)
+    if not numpy.isfinite(da).all() or not numpy.isfinite(db).all() or (da.real <= 0).any() or (db.real <= 0).any():
-    dx_bi = dxes[1][axis].astype(complex)
+        raise ValueError('cell widths must be finite with positive real parts')
-
+    position = numpy.concatenate(([0.0], da.real.cumsum()))
-    if thickness == 0:
+    electric = (position[:-1] + position[1:]) / 2
-        dxes[0][axis] = dx_ai
+    magnetic = position[:-1]
-        dxes[1][axis] = dx_bi
+    if polarity == -1:
-        return dxes
+        selected = slice(0, thickness)
-
+        bound = position[thickness]
-    pos = numpy.hstack((0, dx_ai.cumsum()))
+        length = bound - position[0]
-    pos_a = (pos[:-1] + pos[1:]) / 2
+        depths = [(bound - coordinates[selected]) / length for coordinates in (electric, magnetic)]
-    pos_b = pos[:-1]
-
-    # divide by this to adjust for epsilon_effective and omega
-    s_correction = numpy.sqrt(epsilon_effective) * numpy.real(omega)
-
-    if polarity > 0:
-        # front pml
-        bound = pos[thickness]
-        d = bound - pos[0]
-
-        def l_d(x: NDArray[numpy.float64]) -> NDArray[numpy.float64]:
-            return (bound - x) / (bound - pos[0])
-
-        slc = slice(thickness)
-
     else:
-        # back pml
+        selected = slice(-thickness, None)
-        bound = pos[-thickness - 1]
+        bound = position[-thickness - 1]
-        d = pos[-1] - bound
+        length = position[-1] - bound
-
+        depths = [(coordinates[selected] - bound) / length for coordinates in (electric, magnetic)]
-        def l_d(x: NDArray[numpy.float64]) -> NDArray[numpy.float64]:
+    correction = length * numpy.sqrt(epsilon_effective) * numpy.real(omega)
-            return (x - bound) / (pos[-1] - bound)
+    da[selected] *= _stretch_factors(depths[0], correction, s_function, kappa_function, max_imag_ratio)
-
+    db[selected] *= _stretch_factors(depths[1], correction, s_function, kappa_function, max_imag_ratio)
-        slc = slice(-thickness, None)
+    widths[0][axis], widths[1][axis] = da, db
-
+    return widths
-    dx_ai[slc] *= 1 + 1j * s_function(l_d(pos_a[slc])) / d / s_correction
-    dx_bi[slc] *= 1 + 1j * s_function(l_d(pos_b[slc])) / d / s_correction
-
-    dxes[0][axis] = dx_ai
-    dxes[1][axis] = dx_bi
-
-    return dxes

meanas/fdfd/waveguide_3d.py

@@ -10,6 +10,8 @@ The intended workflow is:
 2. Solve the corresponding 2D mode problem with `solve_mode(...)`.
 3. Turn that mode into a one-sided source with `compute_source(...)`.
 4. Build an overlap window with `compute_overlap_e(...)` for port readout.
+5. For two-port discontinuities, `get_s_fdfd(...)` automates the source,
+   solve, and monitor-plane modal decomposition loop.
 
 `polarity` is part of the public convention throughout this module:
 
@@ -19,15 +21,16 @@ The intended workflow is:
 That same convention controls which side of the selected slice is used for the
 overlap window and how the expanded field is phased.
 """
-from typing import Any, TypedDict, cast
+from typing import Any, TypedDict, cast, Literal
 import warnings
-from collections.abc import Sequence
+import operator
+from collections.abc import Callable, Sequence
 import numpy
 from numpy.typing import NDArray
-from numpy import complexfloating
+from numpy import complexfloating, floating
 
-from ..fdmath import vec, unvec, dx_lists_t, cfdfield_t, fdfield, cfdfield
+from ..fdmath import vec, unvec, dx_lists2_t, dx_lists_t, cfdfield_t, fdfield, cfdfield
-from . import operators, waveguide_2d
+from . import functional, operators, scpml, solvers, waveguide_2d
 
 
 class Waveguide3DMode(TypedDict):
@@ -75,40 +78,86 @@ def solve_mode(
         The returned fields are normalized through the `waveguide_2d`
         normalization convention before being expanded back to 3D.
     """
+    args_2d, _order, _reverse, _dx_prop = _mode_cross_section(omega, dxes, axis, slices, epsilon, mu)
+    e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d)
+    return mode_from_2d(e_xy, wavenumber_2d, omega=omega, dxes=dxes, axis=axis,
+                        polarity=polarity, slices=slices, epsilon=epsilon, mu=mu)
+
+
+def _mode_cross_section(
+        omega: complex,
+        dxes: dx_lists_t,
+        axis: int,
+        slices: Sequence[slice],
+        epsilon: fdfield,
+        mu: fdfield | None,
+        ) -> tuple[dict[str, Any], NDArray, NDArray, complex]:
+    if axis not in range(3):
+        raise ValueError('axis must be 0, 1, or 2')
+    if numpy.ndim(epsilon) != 4 or epsilon.shape[0] != 3:
+        raise ValueError('epsilon must have shape (3, nx, ny, nz)')
+    slices = tuple(slices)
+    if len(slices) != 3 or any(not isinstance(part, slice) for part in slices):
+        raise ValueError('slices must contain three slices')
+    sizes = [len(range(size)[part]) for size, part in zip(epsilon.shape[1:], slices, strict=True)]
+    if sizes[axis] != 1 or any(size == 0 for size in sizes):
+        raise ValueError('slices must select one propagation cell and a nonempty cross-section')
     if mu is None:
         mu = numpy.ones_like(epsilon)
-
+    elif mu.shape != epsilon.shape:
-    slices = tuple(slices)
+        raise ValueError('mu must match epsilon')
-
-    #
-    # Solve the 2D problem in the specified plane
-    #
-    # Define rotation to set z as propagation direction
     order = numpy.roll(range(3), 2 - axis)
-    reverse_order = numpy.roll(range(3), axis - 2)
+    reverse = numpy.roll(range(3), axis - 2)
-
+    dx_prop = 0.5 * sum(dx[axis][slices[axis]].sum() for dx in dxes)
-    # Find dx in propagation direction
+    if not numpy.isfinite(dx_prop) or dx_prop == 0:
-    dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes])
+        raise ValueError('propagation cell width must be finite and nonzero')
-    dx_prop = 0.5 * dxab_forward.sum()
+    args: dict[str, Any] = {
-
-    # Reduce to 2D and solve the 2D problem
-    args_2d: dict[str, Any] = {
         'omega': omega,
         'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes],
         'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]),
         'mu': vec([mu[i][slices].transpose(order) for i in order]),
     }
-    e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d)
+    return args, order, reverse, dx_prop
 
-    #
-    # Apply corrections and expand to 3D
-    #
 
+def mode_from_2d(
+        e_xy: NDArray,
+        wavenumber_2d: complex,
+        *,
+        omega: complex,
+        dxes: dx_lists_t,
+        axis: int,
+        polarity: int,
+        slices: Sequence[slice],
+        epsilon: fdfield,
+        mu: fdfield | None = None,
+        normalization: Literal['power', 'lorentz'] = 'power',
+        ) -> Waveguide3DMode:
+    """Convert a supplied 2D electric eigenmode to a staggered full-grid mode.
+
+    ``e_xy`` contains vectorized Ex/Ey in the rotated plane (the order used by
+    solve_mode). It must solve the 2D problem extracted from ``epsilon`` and
+    ``dxes`` at ``slices``. No eigensolver is called. The helper reconstructs E/H,
+    normalizes them, applies propagation-grid dispersion and half-cell shifts,
+    rotates components, and embeds them in the full epsilon-shaped grid.
+
+    ``polarity`` is +1 toward increasing indices and -1 toward decreasing
+    indices. ``normalization`` follows waveguide_2d.normalized_fields_e.
+    """
+    if polarity not in (-1, 1):
+        raise ValueError('polarity must be -1 or +1')
+    args_2d, _order, reverse_order, dx_prop = _mode_cross_section(omega, dxes, axis, slices, epsilon, mu)
+    slices = tuple(slices)
+    expected = 2 * numpy.prod([len(values) for values in args_2d['dxes'][0]])
+    e_xy = numpy.asarray(e_xy)
+    if e_xy.shape != (expected,) or not numpy.isfinite(e_xy).all() or not numpy.isfinite(wavenumber_2d):
+        raise ValueError('the 2D eigenvector and wavenumber must be finite and match the cross-section')
+    wavenumber_2d = complex(wavenumber_2d)
     # Correct wavenumber to account for numerical dispersion.
     wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2)
 
     shape = [d.size for d in args_2d['dxes'][0]]
-    ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, prop_phase=dx_prop * wavenumber, **args_2d)
+    ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, prop_phase=dx_prop * wavenumber, normalization=normalization, **args_2d)
     e = unvec(ve, shape)
     h = unvec(vh, shape)
 
@@ -160,6 +209,8 @@ def compute_source(
         polarity: Propagation direction (+1 for +ve, -1 for -ve)
         slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
                 as the waveguide cross-section. `slices[axis]` should select only one item.
+        epsilon: Dielectric distribution on the same full grid as E and dxes.
+            Used by the boundary-source operator.
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
@@ -341,3 +392,552 @@ def expand_e(
 
     E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
     return cfdfield_t(E_expanded)
+
+
+def _slice_start(name: str, slices: Sequence[slice], axis: int, shape: Sequence[int]) -> int:
+    if len(slices) != 3:
+        raise ValueError(f'{name} must contain exactly three slices')
+    axis_slice = slices[axis]
+    if axis_slice.start is None or axis_slice.stop is None:
+        raise ValueError(f'{name}[axis] must have explicit start and stop')
+    if axis_slice.step not in (None, 1):
+        raise ValueError(f'{name}[axis] must not use a step')
+    if axis_slice.stop - axis_slice.start != 1:
+        raise ValueError(f'{name}[axis] must select exactly one cell')
+    if axis_slice.start < 0 or axis_slice.stop > shape[axis]:
+        raise ValueError(f'{name}[axis] lies outside the simulation domain')
+    return axis_slice.start
+
+
+def _as_slice_tuple(name: str, slices: Sequence[slice], axis: int, shape: Sequence[int]) -> tuple[slice, slice, slice]:
+    _slice_start(name, slices, axis, shape)
+    return (slices[0], slices[1], slices[2])
+
+
+def _nonbool_index(value: object) -> int:
+    if isinstance(value, bool | numpy.bool_):
+        raise TypeError
+    return operator.index(value)
+
+
+def _as_mode_numbers(name: str, mode_numbers: Sequence[int]) -> tuple[int, ...]:
+    modes: list[int] = []
+    for mode in mode_numbers:
+        try:
+            mode_i = _nonbool_index(mode)
+        except TypeError as err:
+            raise ValueError(f'{name} entries must be integers') from err
+        if mode_i < 0:
+            raise ValueError('mode numbers must be non-negative')
+        modes.append(mode_i)
+    if not modes:
+        raise ValueError('each FDFD port must contain at least one mode number')
+    return tuple(modes)
+
+
+def _validate_fdfd_port_specs(
+        left_mode_numbers: Sequence[int],
+        right_mode_numbers: Sequence[int],
+        left_source_slice: Sequence[slice],
+        left_monitor_slice: Sequence[slice],
+        right_source_slice: Sequence[slice],
+        right_monitor_slice: Sequence[slice],
+        axis: int,
+        shape: Sequence[int],
+        ) -> tuple[
+            tuple[int, ...],
+            tuple[int, ...],
+            tuple[slice, slice, slice],
+            tuple[slice, slice, slice],
+            tuple[slice, slice, slice],
+            tuple[slice, slice, slice],
+            ]:
+    if axis not in (0, 1, 2):
+        raise ValueError('axis must be 0, 1, or 2')
+
+    left_modes_i = _as_mode_numbers('left_mode_numbers', left_mode_numbers)
+    right_modes_i = _as_mode_numbers('right_mode_numbers', right_mode_numbers)
+
+    left_source = _as_slice_tuple('left_source_slice', left_source_slice, axis, shape)
+    left_monitor = _as_slice_tuple('left_monitor_slice', left_monitor_slice, axis, shape)
+    right_source = _as_slice_tuple('right_source_slice', right_source_slice, axis, shape)
+    right_monitor = _as_slice_tuple('right_monitor_slice', right_monitor_slice, axis, shape)
+
+    left_source_i = _slice_start('left_source_slice', left_source, axis, shape)
+    left_monitor_i = _slice_start('left_monitor_slice', left_monitor, axis, shape)
+    right_monitor_i = _slice_start('right_monitor_slice', right_monitor, axis, shape)
+    right_source_i = _slice_start('right_source_slice', right_source, axis, shape)
+    if not (left_source_i < left_monitor_i < right_monitor_i < right_source_i):
+        raise ValueError('port slices must be ordered as left source, left monitor, right monitor, right source')
+
+    return left_modes_i, right_modes_i, left_source, left_monitor, right_source, right_monitor
+
+
+def _port_decomposition(
+        incoming_modes: Sequence[Waveguide3DMode],
+        outgoing_modes: Sequence[Waveguide3DMode],
+        *,
+        dxes: dx_lists_t,
+        axis: int,
+        incoming_polarity: int,
+        outgoing_polarity: int,
+        slices: Sequence[slice],
+        ) -> tuple[NDArray[numpy.complex128], tuple[cfdfield_t, ...], int]:
+    mode_count = len(incoming_modes)
+    if len(outgoing_modes) != mode_count:
+        raise ValueError('incoming and outgoing monitor bases must contain the same number of modes')
+    overlaps = (
+        *(
+            compute_overlap_e(
+                E=mode['E'],
+                wavenumber=mode['wavenumber'],
+                dxes=dxes,
+                axis=axis,
+                polarity=incoming_polarity,
+                slices=slices,
+            )
+            for mode in incoming_modes
+        ),
+        *(
+            compute_overlap_e(
+                E=mode['E'],
+                wavenumber=mode['wavenumber'],
+                dxes=dxes,
+                axis=axis,
+                polarity=outgoing_polarity,
+                slices=slices,
+            )
+            for mode in outgoing_modes
+        ),
+    )
+    basis_fields = (
+        *(
+            expand_e(
+                E=mode['E'],
+                wavenumber=mode['wavenumber'],
+                dxes=dxes,
+                axis=axis,
+                polarity=incoming_polarity,
+                slices=slices,
+            )
+            for mode in incoming_modes
+        ),
+        *(
+            expand_e(
+                E=mode['E'],
+                wavenumber=mode['wavenumber'],
+                dxes=dxes,
+                axis=axis,
+                polarity=outgoing_polarity,
+                slices=slices,
+            )
+            for mode in outgoing_modes
+        ),
+    )
+    matrix = numpy.empty((2 * mode_count, 2 * mode_count), dtype=complex)
+    for row, overlap in enumerate(overlaps):
+        overlap_vec_conj = vec(overlap).conj()
+        for col, basis_field in enumerate(basis_fields):
+            matrix[row, col] = vec(basis_field) @ overlap_vec_conj
+    return matrix, overlaps, mode_count
+
+
+def _decompose_port_field(
+        decomposition: tuple[NDArray[numpy.complex128], tuple[cfdfield_t, ...], int],
+        field: cfdfield_t,
+        ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
+    matrix, overlaps, mode_count = decomposition
+    field_vec = vec(field)
+    measured = numpy.asarray([field_vec @ vec(overlap).conj() for overlap in overlaps])
+    coeffs, _resid, _rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
+    return coeffs[:mode_count], coeffs[mode_count:]
+
+
+def _copy_cross_section(
+        target: NDArray,
+        cross_section: NDArray,
+        axis: int,
+        index: int,
+        ) -> None:
+    slices: list[slice | int] = [slice(None), slice(None), slice(None), slice(None)]
+    slices[axis + 1] = index
+    target[tuple(slices)] = cross_section
+
+
+def _validate_cross_section(name: str, field: NDArray) -> NDArray:
+    array = numpy.asarray(field)
+    if array.ndim != 3 or array.shape[0] != 3:
+        raise ValueError(f'{name} must have shape (3, n0, n1)')
+    if not numpy.isfinite(array).all():
+        raise ValueError(f'{name} must contain only finite values')
+    return array
+
+
+def _as_axis_pair(name: str, value: int | Sequence[int]) -> tuple[int, int]:
+    if isinstance(value, int):
+        pair = (value, value)
+    else:
+        pair = tuple(operator.index(v) for v in value)
+        if len(pair) != 2:
+            raise ValueError(f'{name} must be an integer or a two-entry sequence')
+    if min(pair) < 0:
+        raise ValueError(f'{name} entries must be non-negative')
+    return pair
+
+
+def _pad_cross_section(field: NDArray, padding_cells: tuple[int, int]) -> NDArray:
+    if padding_cells == (0, 0):
+        return field
+    return numpy.pad(
+        field,
+        (
+            (0, 0),
+            (padding_cells[0], padding_cells[0]),
+            (padding_cells[1], padding_cells[1]),
+        ),
+        mode='edge',
+    )
+
+
+def _pad_dxes_2d(dxes_2d: dx_lists2_t, padding_cells: tuple[int, int]) -> list[list[NDArray[floating | complexfloating]]]:
+    padded: list[list[NDArray[floating | complexfloating]]] = [[numpy.empty(0), numpy.empty(0)] for _ in range(2)]
+    for eh in range(2):
+        for axis_2d in range(2):
+            values = numpy.asarray(dxes_2d[eh][axis_2d])
+            pad = padding_cells[axis_2d]
+            if pad == 0:
+                padded[eh][axis_2d] = values.copy()
+            else:
+                padded[eh][axis_2d] = numpy.pad(values, (pad, pad), mode='edge')
+    return padded
+
+
+def get_s_fdfd(
+        omega: complex,
+        dxes: dx_lists_t,
+        epsilon: fdfield,
+        *,
+        left_mode_numbers: Sequence[int],
+        left_source_slice: Sequence[slice],
+        left_monitor_slice: Sequence[slice],
+        right_source_slice: Sequence[slice],
+        right_monitor_slice: Sequence[slice],
+        right_mode_numbers: Sequence[int] | None = None,
+        axis: int = 0,
+        solve_dxes: dx_lists_t | None = None,
+        mu: fdfield | None = None,
+        matrix_solver: Callable[..., Any] | None = None,
+        matrix_solver_opts: dict[str, Any] | None = None,
+        source_scales: Sequence[complex] | None = None,
+        incident_amplitude_atol: float = 1e-12,
+        return_diagnostics: bool = False,
+        ) -> NDArray[numpy.complex128] | tuple[NDArray[numpy.complex128], dict[str, Any]]:
+    """
+    Extract a two-port guided-mode scattering matrix from local 3D FDFD solves.
+
+    The returned S-matrix is ordered as `[[R_L, T_R_to_L], [T_L_to_R, R_R]]`, with
+    reference planes at the supplied monitor slices. Each column is normalized
+    by the incident modal amplitude measured on the corresponding input monitor.
+    """
+    if incident_amplitude_atol <= 0:
+        raise ValueError('incident_amplitude_atol must be positive')
+
+    epsilon = numpy.asarray(epsilon)
+    if epsilon.ndim != 4 or epsilon.shape[0] != 3:
+        raise ValueError('epsilon must have shape (3, nx, ny, nz)')
+    if mu is not None:
+        mu = numpy.asarray(mu)
+        if mu.shape != epsilon.shape:
+            raise ValueError('mu must have the same shape as epsilon')
+
+    if right_mode_numbers is None:
+        right_mode_numbers = left_mode_numbers
+
+    shape = epsilon.shape[1:]
+    left_modes_i, right_modes_i, left_source, left_monitor, right_source, right_monitor = _validate_fdfd_port_specs(
+        left_mode_numbers,
+        right_mode_numbers,
+        left_source_slice,
+        left_monitor_slice,
+        right_source_slice,
+        right_monitor_slice,
+        axis,
+        shape,
+    )
+    n_left = len(left_modes_i)
+    n_right = len(right_modes_i)
+    n_total = n_left + n_right
+
+    if solve_dxes is None:
+        solve_dxes = dxes
+    if source_scales is None:
+        source_scale_array = numpy.ones((n_total,), dtype=complex)
+    else:
+        source_scale_array = numpy.asarray(source_scales, dtype=complex)
+        if source_scale_array.shape != (n_total,):
+            raise ValueError('source_scales must have one entry per incident mode')
+
+    def solve_modes(mode_numbers: Sequence[int], polarity: int, slices: Sequence[slice]) -> tuple[Waveguide3DMode, ...]:
+        return tuple(
+            solve_mode(mode_number, omega=omega, dxes=dxes, axis=axis, polarity=polarity, slices=slices, epsilon=epsilon, mu=mu)
+            for mode_number in mode_numbers
+        )
+
+    left_source_modes = solve_modes(left_modes_i, +1, left_source)
+    right_source_modes = solve_modes(right_modes_i, -1, right_source)
+
+    left_in_modes = solve_modes(left_modes_i, +1, left_monitor)
+    left_out_modes = solve_modes(left_modes_i, -1, left_monitor)
+    right_in_modes = solve_modes(right_modes_i, -1, right_monitor)
+    right_out_modes = solve_modes(right_modes_i, +1, right_monitor)
+    left_decomposition = _port_decomposition(
+        left_in_modes,
+        left_out_modes,
+        dxes=dxes,
+        axis=axis,
+        incoming_polarity=+1,
+        outgoing_polarity=-1,
+        slices=left_monitor,
+    )
+    right_decomposition = _port_decomposition(
+        right_in_modes,
+        right_out_modes,
+        dxes=dxes,
+        axis=axis,
+        incoming_polarity=-1,
+        outgoing_polarity=+1,
+        slices=right_monitor,
+    )
+
+    s = numpy.zeros((n_total, n_total), dtype=complex)
+    raw_outgoing = numpy.zeros_like(s)
+    incident_amplitudes = numpy.zeros((n_total,), dtype=complex)
+    electric_fields: list[cfdfield_t] = []
+    magnetic_fields: list[cfdfield_t] = []
+
+    excitation_modes: list[tuple[int, int, Sequence[slice], Waveguide3DMode]] = [
+        *((col, +1, left_source, left_source_modes[col]) for col in range(n_left)),
+        *((n_left + col, -1, right_source, right_source_modes[col]) for col in range(n_right)),
+    ]
+
+    for column, polarity, source_slice, source_mode in excitation_modes:
+        current = compute_source(
+            E=source_mode['E'],
+            wavenumber=source_mode['wavenumber'],
+            omega=omega,
+            dxes=dxes,
+            axis=axis,
+            polarity=polarity,
+            slices=source_slice,
+            epsilon=epsilon,
+            mu=mu,
+        )
+        current = cfdfield_t(source_scale_array[column] * current)
+        solve_kwargs: dict[str, Any] = {
+            'omega': omega,
+            'dxes': solve_dxes,
+            'J': vec(current),
+            'epsilon': vec(epsilon),
+            'mu': vec(mu),
+            'matrix_solver_opts': {} if matrix_solver_opts is None else dict(matrix_solver_opts),
+        }
+        if matrix_solver is not None:
+            solve_kwargs['matrix_solver'] = matrix_solver
+        e_field = unvec(solvers.generic(**solve_kwargs), shape)
+        left_in, left_out = _decompose_port_field(left_decomposition, e_field)
+        right_in, right_out = _decompose_port_field(right_decomposition, e_field)
+        if column < n_left:
+            incident = left_in[column]
+        else:
+            incident = right_in[column - n_left]
+        if abs(incident) <= incident_amplitude_atol:
+            raise ValueError(f'incident modal amplitude for column {column} is zero')
+
+        raw_column = numpy.concatenate((left_out, right_out))
+
+        incident_amplitudes[column] = incident
+        raw_outgoing[:, column] = raw_column
+        s[:, column] = raw_column / incident
+        if return_diagnostics:
+            electric_fields.append(cfdfield_t(e_field))
+            magnetic_fields.append(cfdfield_t(functional.e2h(omega, solve_dxes, mu=mu)(e_field)))
+
+    if not return_diagnostics:
+        return s
+
+    guided_powers = numpy.sum(numpy.abs(s) ** 2, axis=0).real
+    diagnostics: dict[str, Any] = {
+        'axis': axis,
+        'left_mode_numbers': left_modes_i,
+        'right_mode_numbers': right_modes_i,
+        'left_source_slice': left_source,
+        'left_monitor_slice': left_monitor,
+        'right_source_slice': right_source,
+        'right_monitor_slice': right_monitor,
+        'left_wavenumbers': numpy.asarray([mode['wavenumber'] for mode in left_out_modes]),
+        'right_wavenumbers': numpy.asarray([mode['wavenumber'] for mode in right_out_modes]),
+        'source_scales': source_scale_array,
+        'electric_fields': tuple(electric_fields),
+        'magnetic_fields': tuple(magnetic_fields),
+        'incident_amplitudes': incident_amplitudes,
+        'raw_outgoing_amplitudes': raw_outgoing,
+        'guided_powers': guided_powers,
+        'unresolved_powers': 1.0 - guided_powers,
+        'left_modes': left_out_modes,
+        'right_modes': right_out_modes,
+    }
+    return s, diagnostics
+
+
+def get_s_fdfd_from_2d_ports(
+        left_epsilon: NDArray,
+        right_epsilon: NDArray,
+        *,
+        omega: complex,
+        dxes_2d: dx_lists2_t,
+        mode_numbers: Sequence[int],
+        right_mode_numbers: Sequence[int] | None = None,
+        axis: int = 0,
+        left_mu: NDArray | None = None,
+        right_mu: NDArray | None = None,
+        prop_dx: float | None = None,
+        pml_cells: int = 8,
+        pml_buffer_cells: int = 2,
+        source_monitor_cells: int = 5,
+        monitor_interface_cells: int = 6,
+        transverse_padding_cells: int | Sequence[int] = 0,
+        transverse_pml_cells: int | Sequence[int] = 0,
+        pml_epsilon_effective: float = 1.0,
+        matrix_solver: Callable[..., Any] | None = None,
+        matrix_solver_opts: dict[str, Any] | None = None,
+        source_scales: Sequence[complex] | None = None,
+        return_diagnostics: bool = False,
+        ) -> NDArray[numpy.complex128] | tuple[NDArray[numpy.complex128], dict[str, Any]]:
+    """
+    Extrude two 2D port cross-sections into a short 3D FDFD interface problem.
+    """
+    if axis not in (0, 1, 2):
+        raise ValueError('axis must be 0, 1, or 2')
+    if min(pml_cells, pml_buffer_cells, source_monitor_cells, monitor_interface_cells) < 1:
+        raise ValueError('PML, buffer, source/monitor, and monitor/interface cell counts must be positive')
+
+    padding_cells = _as_axis_pair('transverse_padding_cells', transverse_padding_cells)
+    transverse_pml = _as_axis_pair('transverse_pml_cells', transverse_pml_cells)
+
+    left_epsilon = _pad_cross_section(_validate_cross_section('left_epsilon', left_epsilon), padding_cells)
+    right_epsilon = _pad_cross_section(_validate_cross_section('right_epsilon', right_epsilon), padding_cells)
+    if left_epsilon.shape != right_epsilon.shape:
+        raise ValueError('left_epsilon and right_epsilon must have the same shape')
+    if left_mu is not None:
+        left_mu = _pad_cross_section(_validate_cross_section('left_mu', left_mu), padding_cells)
+        if left_mu.shape != left_epsilon.shape:
+            raise ValueError('left_mu must match left_epsilon')
+    if right_mu is not None:
+        right_mu = _pad_cross_section(_validate_cross_section('right_mu', right_mu), padding_cells)
+        if right_mu.shape != right_epsilon.shape:
+            raise ValueError('right_mu must match right_epsilon')
+    if (left_mu is None) != (right_mu is None):
+        raise ValueError('left_mu and right_mu must be supplied together')
+
+    if right_mode_numbers is None:
+        right_mode_numbers = mode_numbers
+    dxes_2d_padded = _pad_dxes_2d(dxes_2d, padding_cells)
+    transverse_shape = left_epsilon.shape[1:]
+    transverse_axes = tuple(ax for ax in (0, 1, 2) if ax != axis)
+    if prop_dx is None:
+        prop_dx = float(numpy.median(numpy.real(numpy.concatenate([numpy.asarray(d) for d in dxes_2d_padded[0]]))))
+    if prop_dx <= 0:
+        raise ValueError('prop_dx must be positive')
+
+    left_source_i = pml_cells + pml_buffer_cells
+    left_monitor_i = left_source_i + source_monitor_cells
+    interface_i = left_monitor_i + monitor_interface_cells
+    right_monitor_i = interface_i + monitor_interface_cells
+    right_source_i = right_monitor_i + source_monitor_cells
+    axis_size = right_source_i + pml_buffer_cells + pml_cells + 1
+
+    shape3 = [0, 0, 0]
+    shape3[axis] = axis_size
+    shape3[transverse_axes[0]] = transverse_shape[0]
+    shape3[transverse_axes[1]] = transverse_shape[1]
+
+    epsilon = numpy.empty((3, *shape3), dtype=numpy.result_type(left_epsilon, right_epsilon))
+    mu = None if left_mu is None else numpy.empty((3, *shape3), dtype=numpy.result_type(left_mu, right_mu))
+    for index in range(axis_size):
+        cross_epsilon = left_epsilon if index < interface_i else right_epsilon
+        _copy_cross_section(epsilon, cross_epsilon, axis, index)
+        if mu is not None:
+            assert left_mu is not None
+            assert right_mu is not None
+            cross_mu = left_mu if index < interface_i else right_mu
+            _copy_cross_section(mu, cross_mu, axis, index)
+
+    dxes: list[list[NDArray[floating | complexfloating]]] = [[numpy.empty(0), numpy.empty(0), numpy.empty(0)] for _ in range(2)]
+    for eh in range(2):
+        dxes[eh][axis] = numpy.full(axis_size, prop_dx, dtype=float)
+        dxes[eh][transverse_axes[0]] = dxes_2d_padded[eh][0].copy()
+        dxes[eh][transverse_axes[1]] = dxes_2d_padded[eh][1].copy()
+
+    solve_dxes: dx_lists_t = [[d.copy() for d in dx] for dx in dxes]
+    for stretch_axis, thickness in ((axis, pml_cells), (transverse_axes[0], transverse_pml[0]), (transverse_axes[1], transverse_pml[1])):
+        if thickness == 0:
+            continue
+        solve_dxes = scpml.stretch_with_scpml(
+            solve_dxes,
+            axis=stretch_axis,
+            polarity=+1,
+            omega=float(numpy.real(omega)),
+            epsilon_effective=pml_epsilon_effective,
+            thickness=thickness,
+        )
+        solve_dxes = scpml.stretch_with_scpml(
+            solve_dxes,
+            axis=stretch_axis,
+            polarity=-1,
+            omega=float(numpy.real(omega)),
+            epsilon_effective=pml_epsilon_effective,
+            thickness=thickness,
+        )
+
+    def port_slice(index: int) -> tuple[slice, slice, slice]:
+        slices_l = [slice(None), slice(None), slice(None)]
+        slices_l[axis] = slice(index, index + 1)
+        return (slices_l[0], slices_l[1], slices_l[2])
+
+    result = get_s_fdfd(
+        omega=omega,
+        dxes=dxes,
+        solve_dxes=solve_dxes,
+        epsilon=epsilon,
+        mu=mu,
+        axis=axis,
+        left_mode_numbers=mode_numbers,
+        left_source_slice=port_slice(left_source_i),
+        left_monitor_slice=port_slice(left_monitor_i),
+        right_mode_numbers=right_mode_numbers,
+        right_source_slice=port_slice(right_source_i),
+        right_monitor_slice=port_slice(right_monitor_i),
+        matrix_solver=matrix_solver,
+        matrix_solver_opts=matrix_solver_opts,
+        source_scales=source_scales,
+        return_diagnostics=return_diagnostics,
+    )
+    if not return_diagnostics:
+        return result
+
+    s, diagnostics = result
+    diagnostics.update({
+        'interface_index': interface_i,
+        'pml_cells': pml_cells,
+        'pml_buffer_cells': pml_buffer_cells,
+        'source_monitor_cells': source_monitor_cells,
+        'monitor_interface_cells': monitor_interface_cells,
+        'transverse_padding_cells': padding_cells,
+        'transverse_pml_cells': transverse_pml,
+        'prop_dx': prop_dx,
+        'epsilon': epsilon,
+        'mu': mu,
+        'dxes': dxes,
+        'solve_dxes': solve_dxes,
+    })
+    return s, diagnostics

meanas/fdmath/functional.py

@@ -13,52 +13,61 @@ from numpy import floating, complexfloating
 from .types import fdfield, fdfield_updater_t
 
 
+def _derivative(axis: int, forward: bool, widths: NDArray | None) -> fdfield_updater_t:
+    """Periodic difference using views, with one output allocation at most."""
+    def derivative(f: NDArray, *, out: NDArray | None = None) -> NDArray:
+        f = numpy.asarray(f)
+        if f.ndim != 3 or f.shape[axis] == 0:
+            raise ValueError('derivatives require a nonempty three-dimensional axis')
+        if widths is not None and (numpy.ndim(widths) != 1 or len(widths) != f.shape[axis]):
+            raise ValueError('cell widths must match the differentiated axis')
+        dtype = f.dtype if widths is None else numpy.true_divide.resolve_dtypes((f.dtype, numpy.asarray(widths).dtype, None))[2]
+        if out is None:
+            output = numpy.empty(f.shape, dtype=dtype)
+        else:
+            output = numpy.asarray(out)
+            if output.shape != f.shape:
+                raise ValueError('out must match the field shape')
+            if numpy.shares_memory(f, output) or (widths is not None and numpy.shares_memory(widths, output)):
+                raise ValueError('out must not alias the field or widths')
+            if not output.flags.writeable or not numpy.can_cast(dtype, output.dtype, casting='same_kind'):
+                raise ValueError('out must be writable with a compatible dtype')
+        low, high = [slice(None)] * 3, [slice(None)] * 3
+        low[axis], high[axis] = slice(None, -1), slice(1, None)
+        first, last = [slice(None)] * 3, [slice(None)] * 3
+        first[axis], last[axis] = slice(0, 1), slice(-1, None)
+        destination = low if forward else high
+        numpy.subtract(f[tuple(high)], f[tuple(low)], out=output[tuple(destination)])
+        destination = last if forward else first
+        numpy.subtract(f[tuple(first)], f[tuple(last)], out=output[tuple(destination)])
+        if widths is not None:
+            expand = [1, 1, 1]
+            expand[axis] = len(widths)
+            numpy.divide(output, numpy.asarray(widths).reshape(expand), out=output)
+        return output
+    return derivative
+
+
 def deriv_forward(
         dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]:
-    """
+    """Forward periodic derivatives divided by electric Yee widths.
-    Utility operators for taking discretized derivatives (backward variant).
 
-    Args:
+    Each callable accepts ``derivative(field, *, out=None)``. An output buffer
-        dx_e: Lists of cell sizes for all axes
+    must match the field shape and must not alias the field or width arrays.
-              `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+    Slice-based subtraction avoids a full-size roll temporary.
-
-    Returns:
-        List of functions for taking forward derivatives along each axis.
     """
-    if dx_e is not None:
+    return tuple(_derivative(axis, True, None if dx_e is None else dx_e[axis]) for axis in range(3))
-        derivs = (lambda f: (numpy.roll(f, -1, axis=0) - f) / dx_e[0][:, None, None],
-                  lambda f: (numpy.roll(f, -1, axis=1) - f) / dx_e[1][None, :, None],
-                  lambda f: (numpy.roll(f, -1, axis=2) - f) / dx_e[2][None, None, :])
-    else:
-        derivs = (lambda f: numpy.roll(f, -1, axis=0) - f,
-                  lambda f: numpy.roll(f, -1, axis=1) - f,
-                  lambda f: numpy.roll(f, -1, axis=2) - f)
-    return derivs
 
 
 def deriv_back(
         dx_h: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]:
-    """
+    """Backward periodic derivatives divided by magnetic Yee widths.
-    Utility operators for taking discretized derivatives (forward variant).
 
-    Args:
+    Callables accept the same non-aliasing ``out`` buffer as deriv_forward.
-        dx_h: Lists of cell sizes for all axes
-              `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
-
-    Returns:
-        List of functions for taking forward derivatives along each axis.
     """
-    if dx_h is not None:
+    return tuple(_derivative(axis, False, None if dx_h is None else dx_h[axis]) for axis in range(3))
-        derivs = (lambda f: (f - numpy.roll(f, 1, axis=0)) / dx_h[0][:, None, None],
-                  lambda f: (f - numpy.roll(f, 1, axis=1)) / dx_h[1][None, :, None],
-                  lambda f: (f - numpy.roll(f, 1, axis=2)) / dx_h[2][None, None, :])
-    else:
-        derivs = (lambda f: f - numpy.roll(f, 1, axis=0),
-                  lambda f: f - numpy.roll(f, 1, axis=1),
-                  lambda f: f - numpy.roll(f, 1, axis=2))
-    return derivs
 
 
 TT = TypeVar('TT', bound='NDArray[floating | complexfloating]')

meanas/fdtd/pml.py

@@ -41,6 +41,7 @@ def cpml_params(
         m: float = 3.5,
         ma: float = 1,
         cfs_alpha: float = 0,
+        dx: float = 1.0,
         ) -> dict[str, Any]:
     """
     Construct the parameter block for one CPML face.
@@ -51,12 +52,17 @@ def cpml_params(
             `+1` for the high-index face).
         dt: Timestep used by the Yee update.
         thickness: Number of Yee cells occupied by the CPML region.
-        ln_R_per_layer: Logarithmic attenuation target per layer.
+        ln_R_per_layer: Finite, strictly negative logarithmic attenuation target
+            per layer. Zero attenuation is not supported; disable a face with
+            ``None`` in the updates_with_cpml parameter table instead.
         epsilon_eff: Effective permittivity used when choosing the CPML scaling.
         mu_eff: Effective permeability used when choosing the CPML scaling.
         m: Polynomial grading exponent for `sigma` and `kappa`.
         ma: Polynomial grading exponent for the complex-frequency shift `alpha`.
         cfs_alpha: Maximum complex-frequency shift parameter.
+        dx: Physical cell width normal to this face, uniform within the CPML.
+            The default assumes unit cells. Arbitrary nonuniform grading is
+            not supported by this parameter constructor.
 
     Returns:
         Dictionary with:
@@ -78,7 +84,13 @@ def cpml_params(
     if epsilon_eff <= 0:
         raise ValueError('epsilon_eff must be positive')
 
-    sigma_max = -ln_R_per_layer / 2 * (m + 1)
+    if not numpy.isfinite(ln_R_per_layer) or ln_R_per_layer >= 0:
+        raise ValueError('ln_R_per_layer must be finite and negative; use None to disable a CPML face')
+
+    if not numpy.isfinite(dx) or dx <= 0:
+        raise ValueError('dx must be finite and positive')
+
+    sigma_max = -ln_R_per_layer / (2 * dx) * (m + 1)
     kappa_max = numpy.sqrt(epsilon_eff * mu_eff)
     alpha_max = cfs_alpha
 
@@ -152,7 +164,7 @@ def updates_with_cpml(
         The closures retain the CPML auxiliary state internally.
     """
 
-    Dfx, Dfy, Dfz = deriv_forward(dxes[1])
+    Dfx, Dfy, Dfz = deriv_forward(dxes[0])
     Dbx, Dby, Dbz = deriv_back(dxes[1])
 
     psi_E: list[list[tuple[Any, Any]]] = [[(None, None) for _ in range(2)] for _ in range(3)]

meanas/test/test_eigensolvers_numerics.py

@@ -1,10 +1,11 @@
 import numpy
 from numpy.linalg import norm
 import pytest
+from scipy import sparse
 
 from ._solver_cases import diagonal_eigen_case
 from .utils import assert_close
-from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve
+from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, scipy_eigs_near, signed_eigensolve
 
 
 def test_rayleigh_quotient_iteration_with_linear_operator() -> None:
@@ -99,3 +100,27 @@ def test_signed_eigensolve_returns_largest_positive_modes() -> None:
 
     assert_close(eigvals, [3.0, 5.0], atol=1e-6)
     assert eigvecs.shape == (4, 2)
+
+
+def test_scipy_eigs_near_returns_paired_bilinear_left_vectors() -> None:
+    operator = sparse.csc_array(numpy.diag([1.0 + 0.1j, 2.0 - 0.2j, 4.0 + 0.3j, 7.0]))
+
+    solution = scipy_eigs_near(operator, target=3.1, how_many=2)
+
+    assert_close(numpy.sort_complex(solution.eigenvalues), numpy.sort_complex([2.0 - 0.2j, 4.0 + 0.3j]))
+    gram = solution.left_eigenvectors.T @ solution.right_eigenvectors
+    assert numpy.linalg.cond(gram) < 2
+    for index, eigenvalue in enumerate(solution.eigenvalues):
+        right = solution.right_eigenvectors[:, index]
+        left = solution.left_eigenvectors[:, index]
+        assert norm(operator @ right - eigenvalue * right) < 1e-10
+        assert norm(operator.T @ left - eigenvalue * left) < 1e-10
+
+
+@pytest.mark.parametrize('parameter', ['tolerance', 'pairing_rtol'])
+@pytest.mark.parametrize('invalid', [numpy.nan, numpy.inf, -numpy.inf])
+def test_scipy_eigs_near_rejects_nonfinite_tolerances(parameter: str, invalid: float) -> None:
+    operator = sparse.csc_array(numpy.diag([1.0, 2.0, 3.0, 4.0]))
+
+    with pytest.raises(ValueError, match=rf'{parameter} must be finite and positive'):
+        scipy_eigs_near(operator, target=2.1, how_many=1, **{parameter: invalid})

meanas/test/test_eme_numerics.py

@@ -2,10 +2,9 @@ from typing import cast
 
 import numpy
 import pytest
-from scipy import sparse
 
 from ..fdmath import vec
-from ..fdfd import eme, waveguide_2d, waveguide_cyl
+from ..fdfd import eme, waveguide_2d, waveguide_3d, waveguide_cyl
 from ._test_builders import complex_ramp, unit_dxes
 from .utils import assert_close
 
@@ -15,7 +14,10 @@ DXES = unit_dxes((2, 2))
 WAVENUMBERS_L = numpy.array([1.0, 0.8])
 WAVENUMBERS_R = numpy.array([0.9, 0.7])
 OMEGA = 1 / 1500
+FDFD_OMEGA = 2 * numpy.pi / 16
 REAL_DXES = unit_dxes((5, 5))
+FDFD_SHAPE_2D = (7, 7)
+FDFD_DXES_2D = unit_dxes(FDFD_SHAPE_2D)
 
 
 def _mode(scale: float) -> tuple[numpy.ndarray, numpy.ndarray]:
@@ -30,34 +32,11 @@ def _mode_sets() -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], list[tuple[
     return left_modes, right_modes
 
 
-def _gain_only_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
-    return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.zeros((2, 2))
-
-
-def _gain_and_reflection_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
-    return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.array([[0.0, 1.0], [2.0, 0.0]])
-
-
-def _nonsymmetric_tr(left_marker: object):
-    def fake_get_tr(_eh_left, wavenumbers_left, _eh_right, _wavenumbers_right, **kwargs):
-        if wavenumbers_left is left_marker:
-            return (
-                numpy.array([[1.0, 2.0], [0.5, 1.0]]),
-                numpy.array([[0.0, 1.0], [2.0, 0.0]]),
-            )
-        return (
-            numpy.array([[1.0, -1.0], [0.0, 1.0]]),
-            numpy.array([[0.0, 0.5], [1.5, 0.0]]),
-        )
-
-    return fake_get_tr
-
-
 def _dummy_modes() -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]:
     return [_mode(0.0), _mode(0.7)], numpy.array([1.0, 0.5])
 
 
-def test_get_tr_returns_finite_bounded_transfer_matrices() -> None:
+def test_get_tr_returns_finite_coordinate_transfer_matrices() -> None:
     left_modes, right_modes = _mode_sets()
 
     transmission, reflection = eme.get_tr(
@@ -74,7 +53,8 @@ def test_get_tr_returns_finite_bounded_transfer_matrices() -> None:
     assert reflection.shape == (2, 2)
     assert numpy.isfinite(transmission).all()
     assert numpy.isfinite(reflection).all()
-    assert (singular_values <= 1.0 + 1e-12).all()
+    # These arbitrary fields are not power-normalized: coefficients may exceed one.
+    assert singular_values.max() > 1.0
 
 
 def test_get_tr_accepts_scaled_dual_projection_modes() -> None:
@@ -95,73 +75,53 @@ def test_get_tr_accepts_scaled_dual_projection_modes() -> None:
     )
 
     assert_close(dual_t, plain_t)
-    assert_close(dual_r, plain_r)
+    assert_close(dual_r, plain_r, atol=1e-11)
 
 
-def test_get_abcd_matches_explicit_block_formula() -> None:
+def test_get_abcd_maps_left_state_to_right_state() -> None:
     left_modes, right_modes = _mode_sets()
-    t12, r12 = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+    ss = eme.get_s(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, DXES)
-    t21, r21 = eme.get_tr(right_modes, WAVENUMBERS_R, left_modes, WAVENUMBERS_L, dxes=DXES)
+    transfer = eme.get_abcd(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, DXES).toarray()
-    t21_inv = numpy.linalg.pinv(t21)
+    incoming = numpy.array([1.0, .3j, -.4j, .2])
-
+    outgoing = ss @ incoming
-    expected = numpy.block([
+    state_left = numpy.concatenate((incoming[:2], outgoing[:2]))
-        [t12 - r21 @ t21_inv @ r12, r21 @ t21_inv],
+    state_right = numpy.concatenate((outgoing[2:], incoming[2:]))
-        [-t21_inv @ r12, t21_inv],
+    assert_close(transfer @ state_left, state_right, atol=1e-11)
-    ])
-    abcd = eme.get_abcd(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
-
-    assert sparse.issparse(abcd)
-    assert abcd.shape == (4, 4)
-    assert_close(abcd.toarray(), expected)
 
 
-def test_get_s_plain_matches_block_assembly_from_get_tr() -> None:
+def test_get_s_left_incidence_matches_get_tr() -> None:
     left_modes, right_modes = _mode_sets()
-    t12, r12 = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+    transmission, reflection = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, DXES)
-    t21, r21 = eme.get_tr(right_modes, WAVENUMBERS_R, left_modes, WAVENUMBERS_L, dxes=DXES)
+    ss = eme.get_s(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, DXES)
-    expected = numpy.block([[r12, t12], [t21, r21]])
-
-    ss = eme.get_s(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
-
     assert ss.shape == (4, 4)
-    assert numpy.isfinite(ss).all()
+    assert_close(ss[:2, :2], reflection, atol=1e-11)
-    assert_close(ss, expected)
+    assert_close(ss[2:, :2], transmission, atol=1e-11)
 
 
 def test_get_s_force_nogain_caps_singular_values(monkeypatch) -> None:
-    monkeypatch.setattr(eme, 'get_tr', _gain_only_tr)
+    raw = numpy.diag([2., .5, 1.2, .3])
+    monkeypatch.setattr(eme, '_interface_s', lambda *args, **kwargs: raw.copy())
     modes, wavenumbers = _dummy_modes()
-
+    plain_s = eme.get_s(modes, wavenumbers, modes, wavenumbers, DXES)
-    plain_s = eme.get_s(modes, wavenumbers, modes, wavenumbers)
+    clipped_s = eme.get_s(modes, wavenumbers, modes, wavenumbers, DXES, force_nogain=True)
-    clipped_s = eme.get_s(modes, wavenumbers, modes, wavenumbers, force_nogain=True)
+    assert numpy.linalg.svd(plain_s, compute_uv=False).max() > 1
-
+    assert_close(numpy.linalg.svd(clipped_s, compute_uv=False), [1., 1., .5, .3])
-    plain_singular_values = numpy.linalg.svd(plain_s, compute_uv=False)
-    clipped_singular_values = numpy.linalg.svd(clipped_s, compute_uv=False)
-
-    assert plain_singular_values.max() > 1.0
-    assert (clipped_singular_values <= 1.0 + 1e-12).all()
-    assert numpy.isfinite(clipped_s).all()
 
 
 def test_get_s_force_reciprocal_symmetrizes_output(monkeypatch) -> None:
-    left = numpy.array([1.0, 0.5])
+    raw = numpy.arange(16, dtype=float).reshape(4, 4)
-    right = numpy.array([0.9, 0.4])
+    monkeypatch.setattr(eme, '_interface_s', lambda *args, **kwargs: raw.copy())
-    modes, _wavenumbers = _dummy_modes()
+    modes, wavenumbers = _dummy_modes()
-
+    ss = eme.get_s(modes, wavenumbers, modes, wavenumbers, DXES, force_reciprocal=True)
-    monkeypatch.setattr(eme, 'get_tr', _nonsymmetric_tr(left))
+    assert_close(ss, .5 * (raw + raw.T))
-    ss = eme.get_s(modes, left, modes, right, force_reciprocal=True)
-
-    assert_close(ss, ss.T)
 
 
 def test_get_s_force_nogain_and_reciprocal_returns_finite_output(monkeypatch) -> None:
-    monkeypatch.setattr(eme, 'get_tr', _gain_and_reflection_tr)
+    raw = numpy.arange(16, dtype=float).reshape(4, 4)
+    monkeypatch.setattr(eme, '_interface_s', lambda *args, **kwargs: raw.copy())
     modes, wavenumbers = _dummy_modes()
-    ss = eme.get_s(modes, wavenumbers, modes, wavenumbers, force_nogain=True, force_reciprocal=True)
+    ss = eme.get_s(modes, wavenumbers, modes, wavenumbers, DXES, force_nogain=True, force_reciprocal=True)
-
+    assert_close(ss, ss.T, atol=1e-14)
-    assert ss.shape == (4, 4)
-    assert numpy.isfinite(ss).all()
-    assert_close(ss, ss.T)
     assert (numpy.linalg.svd(ss, compute_uv=False) <= 1.0 + 1e-12).all()
 
 
@@ -461,8 +421,8 @@ def _abcd_to_s(abcd: numpy.ndarray) -> numpy.ndarray:
     t12 = aa - bb * cc / dd
     return numpy.array(
         [
-            [r12, t12],
+            [r12, t21],
-            [t21, r21],
+            [t12, r21],
         ],
         dtype=complex,
     )
@@ -481,6 +441,25 @@ def _expected_bragg_reflector_s(n_low: float, n_high: float, periods: int) -> nu
     )
 
 
+def _fdfd_cross_section(core_epsilon: float) -> numpy.ndarray:
+    epsilon = numpy.ones((3, *FDFD_SHAPE_2D), dtype=float)
+    epsilon[:, 2:5, 2:5] = core_epsilon
+    return epsilon
+
+
+def _small_fdfd_kwargs() -> dict:
+    return {
+        'omega': FDFD_OMEGA,
+        'dxes_2d': FDFD_DXES_2D,
+        'mode_numbers': [0],
+        'pml_cells': 4,
+        'pml_buffer_cells': 2,
+        'source_monitor_cells': 5,
+        'monitor_interface_cells': 5,
+        'matrix_solver_opts': {'atol': 1e-10, 'rtol': 1e-8},
+    }
+
+
 def test_get_s_is_near_identity_for_identical_solved_straight_modes() -> None:
     mode, wavenumber, _epsilon = _build_straight_mode()
 
@@ -555,6 +534,189 @@ def test_get_s_accepts_complex_leaky_wavenumbers_for_abrupt_interface() -> None:
     assert_close(ss, numpy.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex), atol=1e-12, rtol=1e-12)
 
 
+def test_get_s_fdfd_from_2d_ports_identical_guides_are_nearly_transmissive() -> None:
+    epsilon = _fdfd_cross_section(2.25)
+
+    ss = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, **_small_fdfd_kwargs())
+
+    assert ss.shape == (2, 2)
+    assert numpy.isfinite(ss).all()
+    assert abs(ss[0, 0]) < 0.03
+    assert abs(ss[1, 1]) < 0.04
+    assert abs(abs(ss[1, 0]) - 1.0) < 0.02
+    assert abs(abs(ss[0, 1]) - 1.0) < 0.03
+
+
+def test_get_s_fdfd_from_2d_ports_mismatch_returns_finite_diagnostics() -> None:
+    left = _fdfd_cross_section(2.25)
+    right = _fdfd_cross_section(1.5)
+
+    ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **_small_fdfd_kwargs())
+
+    assert ss.shape == (2, 2)
+    assert numpy.isfinite(ss).all()
+    assert len(diagnostics['electric_fields']) == 2
+    assert len(diagnostics['magnetic_fields']) == 2
+    assert diagnostics['electric_fields'][0].shape == (3, 33, *FDFD_SHAPE_2D)
+    assert numpy.isfinite(diagnostics['electric_fields'][0]).all()
+    assert numpy.isfinite(diagnostics['magnetic_fields'][0]).all()
+    assert numpy.isfinite(diagnostics['guided_powers']).all()
+    assert numpy.isfinite(diagnostics['unresolved_powers']).all()
+    assert numpy.isfinite(diagnostics['raw_outgoing_amplitudes']).all()
+    assert numpy.isfinite(diagnostics['incident_amplitudes']).all()
+    assert (numpy.abs(diagnostics['incident_amplitudes']) > 1e-3).all()
+    assert_close(diagnostics['unresolved_powers'], 1.0 - diagnostics['guided_powers'])
+    assert diagnostics['guided_powers'][0] < 1.05
+
+
+def test_get_s_fdfd_normalization_is_source_scale_invariant() -> None:
+    left = _fdfd_cross_section(2.25)
+    right = _fdfd_cross_section(1.5)
+    kwargs = _small_fdfd_kwargs()
+
+    base, base_diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **kwargs)
+    scaled, scaled_diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, source_scales=[2.0, 0.3], return_diagnostics=True, **kwargs)
+
+    assert_close(numpy.abs(scaled), numpy.abs(base), atol=1e-6, rtol=1e-6)
+    assert abs(abs(scaled_diagnostics['incident_amplitudes'][0] / base_diagnostics['incident_amplitudes'][0]) - 2.0) < 1e-6
+    assert abs(abs(scaled_diagnostics['incident_amplitudes'][1] / base_diagnostics['incident_amplitudes'][1]) - 0.3) < 1e-6
+
+
+def test_get_s_fdfd_from_2d_ports_diagnostics_order_generated_slices() -> None:
+    epsilon = _fdfd_cross_section(2.25)
+
+    _ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, return_diagnostics=True, **_small_fdfd_kwargs())
+
+    left_source = diagnostics['left_source_slice'][0].start
+    left_monitor = diagnostics['left_monitor_slice'][0].start
+    interface_index = diagnostics['interface_index']
+    right_monitor = diagnostics['right_monitor_slice'][0].start
+    right_source = diagnostics['right_source_slice'][0].start
+
+    assert left_source == 6
+    assert left_monitor == 11
+    assert interface_index == 16
+    assert right_monitor == 21
+    assert right_source == 26
+    assert left_source < left_monitor < interface_index < right_monitor < right_source
+
+
+def test_get_s_fdfd_from_2d_ports_transverse_padding_and_pml_updates_diagnostics() -> None:
+    epsilon = _fdfd_cross_section(2.25)
+
+    ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(
+        epsilon,
+        epsilon,
+        transverse_padding_cells=(2, 3),
+        transverse_pml_cells=(1, 2),
+        return_diagnostics=True,
+        **_small_fdfd_kwargs(),
+    )
+
+    assert ss.shape == (2, 2)
+    assert diagnostics['epsilon'].shape == (3, 33, 11, 13)
+    assert diagnostics['transverse_padding_cells'] == (2, 3)
+    assert diagnostics['transverse_pml_cells'] == (1, 2)
+    assert numpy.iscomplexobj(diagnostics['solve_dxes'][0][1])
+    assert numpy.iscomplexobj(diagnostics['solve_dxes'][0][2])
+    assert numpy.abs(numpy.imag(diagnostics['solve_dxes'][0][1][0])) > 0
+    assert numpy.abs(numpy.imag(diagnostics['solve_dxes'][0][2][0])) > 0
+
+
+def test_get_s_fdfd_from_2d_ports_supports_nondefault_axis() -> None:
+    epsilon = _fdfd_cross_section(2.25)
+
+    ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, axis=1, return_diagnostics=True, **_small_fdfd_kwargs())
+
+    assert ss.shape == (2, 2)
+    assert diagnostics['epsilon'].shape == (3, FDFD_SHAPE_2D[0], 33, FDFD_SHAPE_2D[1])
+    assert diagnostics['left_source_slice'][1].start == 6
+    assert numpy.isfinite(ss).all()
+
+
+def test_get_s_fdfd_rejects_bad_port_slice_order() -> None:
+    epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
+    dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
+
+    with pytest.raises(ValueError, match='ordered'):
+        waveguide_3d.get_s_fdfd(
+            FDFD_OMEGA,
+            dxes,
+            epsilon,
+            left_mode_numbers=[0],
+            left_source_slice=(slice(5, 6), slice(None), slice(None)),
+            left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
+            right_source_slice=(slice(10, 11), slice(None), slice(None)),
+            right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
+        )
+
+
+def test_get_s_fdfd_rejects_bad_source_scale_count_before_solving() -> None:
+    epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
+    dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
+
+    with pytest.raises(ValueError, match='source_scales'):
+        waveguide_3d.get_s_fdfd(
+            FDFD_OMEGA,
+            dxes,
+            epsilon,
+            left_mode_numbers=[0],
+            left_source_slice=(slice(2, 3), slice(None), slice(None)),
+            left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
+            right_mode_numbers=[0, 1],
+            right_source_slice=(slice(10, 11), slice(None), slice(None)),
+            right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
+            source_scales=[1.0, 1.0],
+        )
+
+
+def test_get_s_fdfd_from_2d_ports_rejects_invalid_generated_geometry() -> None:
+    epsilon = _fdfd_cross_section(2.25)
+
+    with pytest.raises(ValueError, match='must be positive'):
+        waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, **{**_small_fdfd_kwargs(), 'pml_cells': 0})
+
+
+def test_get_s_fdfd_rejects_noninteger_mode_numbers_before_solving() -> None:
+    epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
+    dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
+
+    with pytest.raises(ValueError, match='mode_numbers'):
+        waveguide_3d.get_s_fdfd(
+            FDFD_OMEGA,
+            dxes,
+            epsilon,
+            left_mode_numbers=[0.5],  # type: ignore[list-item]
+            left_source_slice=(slice(2, 3), slice(None), slice(None)),
+            left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
+            right_source_slice=(slice(10, 11), slice(None), slice(None)),
+            right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
+        )
+
+
+def test_get_s_fdfd_from_2d_ports_preserves_complex_generated_inputs(monkeypatch: pytest.MonkeyPatch) -> None:
+    left = _fdfd_cross_section(2.25)
+    right = _fdfd_cross_section(1.5).astype(complex)
+    right += 0.1j
+    dxes_2d = [
+        [FDFD_DXES_2D[0][0].astype(complex) + 0.2j, FDFD_DXES_2D[0][1].astype(complex) + 0.3j],
+        [FDFD_DXES_2D[1][0].astype(complex) + 0.4j, FDFD_DXES_2D[1][1].astype(complex) + 0.5j],
+    ]
+
+    def fake_get_s_fdfd(*args, **kwargs) -> tuple[numpy.ndarray, dict]:
+        return numpy.eye(2, dtype=complex), {}
+
+    monkeypatch.setattr(waveguide_3d, 'get_s_fdfd', fake_get_s_fdfd)
+
+    kwargs = {**_small_fdfd_kwargs(), 'dxes_2d': dxes_2d}
+    _ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **kwargs)
+
+    assert numpy.iscomplexobj(diagnostics['epsilon'])
+    assert diagnostics['epsilon'][0, diagnostics['interface_index'], 0, 0].imag == pytest.approx(0.1)
+    assert numpy.iscomplexobj(diagnostics['dxes'][0][1])
+    assert diagnostics['dxes'][0][1][0].imag == pytest.approx(0.2)
+
+
 def test_quarter_wave_matching_layer_is_nearly_reflectionless_at_design_frequency() -> None:
     """
     A single quarter-wave matching layer with
@@ -572,7 +734,7 @@ def test_quarter_wave_matching_layer_is_nearly_reflectionless_at_design_frequenc
     _mode_1, beta_1 = _build_uniform_mode(float(n1))
     quarter_wave_length = numpy.pi / (2 * numpy.real(beta_1))
 
-    stack_abcd = interface_01 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_12
+    stack_abcd = interface_12 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_01
     ss = _abcd_to_s(stack_abcd)
 
     assert ss.shape == (2, 2)
@@ -605,7 +767,7 @@ def test_quarter_wave_ar_layer_reduces_reflection_relative_to_abrupt_interface()
     interface_12 = _interface_abcd(n1, n2)
     _mode_1, beta_1 = _build_uniform_mode(float(n1))
     quarter_wave_length = numpy.pi / (2 * numpy.real(beta_1))
-    ar_stack = _abcd_to_s(interface_01 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_12)
+    ar_stack = _abcd_to_s(interface_12 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_01)
 
     abrupt_reflection = abs(abrupt[0, 0])
     abrupt_transmission = abs(abrupt[1, 0])
@@ -636,7 +798,7 @@ def test_half_wave_uniform_slab_restores_unit_transmission_between_matched_media
     _mode_1, beta_1 = _build_uniform_mode(n1)
     half_wave_length = numpy.pi / numpy.real(beta_1)
 
-    stack_abcd = interface_01 @ _propagation_abcd(beta_1, half_wave_length) @ interface_10
+    stack_abcd = interface_10 @ _propagation_abcd(beta_1, half_wave_length) @ interface_01
     ss = _abcd_to_s(stack_abcd)
 
     assert ss.shape == (2, 2)
@@ -694,8 +856,7 @@ def test_quarter_wave_bragg_reflector_matches_closed_form_stopband_response() ->
 
     stack_abcd = numpy.eye(2, dtype=complex)
     for _ in range(periods):
-        stack_abcd = stack_abcd @ interface_lh @ _propagation_abcd(beta_h, quarter_wave_high)
+        stack_abcd = _propagation_abcd(beta_l, quarter_wave_low) @ interface_hl @ _propagation_abcd(beta_h, quarter_wave_high) @ interface_lh @ stack_abcd
-        stack_abcd = stack_abcd @ interface_hl @ _propagation_abcd(beta_l, quarter_wave_low)
     ss = _abcd_to_s(stack_abcd)
     expected = _expected_bragg_reflector_s(n_low, n_high, periods)
 

meanas/test/test_eme_orientation.py

@@ -0,0 +1,210 @@
+"""Characterize truncation bias against exact reciprocal Maxwell interfaces.
+
+These tests intentionally document mixed/Galerkin projection limitations; they
+are not passivity-enforcement tests. Both complete modal bases and an independent
+driven FDFD solve establish the reference for the truncated polarization case.
+"""
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+from scipy.sparse.linalg import spsolve
+
+from ..fdmath import vec, unvec
+from ..fdfd import eme, scpml, solvers, waveguide_2d, waveguide_3d
+
+
+DXES = [[np.ones(1), np.ones(1)] for _ in range(2)]
+
+
+def _rotation(angle):
+    return np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
+
+
+def _polarizations(angle, scales=(1., 1.)):
+    """Two power-orthonormal signed-beta=1 modes on a one-cell transverse grid.
+
+    E_xy = sqrt(2) diag(scales) Q, and J^-1 H_xy = sqrt(2) diag(1/scales) Q.
+    They solve Maxwell with diagonal epsilon=(1/sx^2,1/sy^2,1) and
+    mu=(sy^2,sx^2,1) at omega=1, including nonparallel E/H traces.
+    """
+    rotation = _rotation(angle)
+    electric = np.diag(scales) @ rotation
+    magnetic = np.diag(1 / np.asarray(scales)) @ rotation
+    return [
+        (np.sqrt(2) * np.array([electric[0, i], electric[1, i], 0], complex),
+         np.sqrt(2) * np.array([-magnetic[1, i], magnetic[0, i], 0], complex))
+        for i in range(2)
+    ]
+
+
+def _scattering(left, right, formulation='mixed'):
+    return eme.get_s(left, np.ones(len(left)), right, np.ones(len(right)), DXES, formulation=formulation)
+
+
+def _reorder_ports(matrix):
+    n = matrix.shape[0] // 2
+    order = [*range(n, 2*n), *range(n)]
+    return matrix[np.ix_(order, order)]
+
+
+def _fresnel_reference(left_angle, right_angle, left_scales=(1., 1.), right_scales=(1., 1.)):
+    """Independent exact two-polarization interface in the supplied port bases."""
+    yl, yr = 1 / np.asarray(left_scales)**2, 1 / np.asarray(right_scales)**2
+    reflection = np.diag((yl - yr) / (yl + yr))
+    transmission = np.diag(2 * np.sqrt(yl * yr) / (yl + yr))
+    ql, qr = _rotation(left_angle), _rotation(right_angle)
+    return np.block([[ql.T @ reflection @ ql, ql.T @ transmission @ qr],
+                     [qr.T @ transmission @ ql, -qr.T @ reflection @ qr]])
+
+
+def _trace_least_squares(left, right):
+    """Research comparator: unweighted Euclidean tangential E/H least squares.
+
+    This is deliberately not a public formulation: correcting the isotropic
+    polarization example does not establish reciprocity or passivity generally.
+    """
+    el = np.column_stack([e[:2] for e, _h in left])
+    hl = np.column_stack([h[:2] for _e, h in left])
+    er = np.column_stack([e[:2] for e, _h in right])
+    hr = np.column_stack([h[:2] for _e, h in right])
+    outgoing = np.block([[el, -er], [-hl, -hr]])
+    incoming = np.block([[-el, er], [-hl, -hr]])
+    return np.linalg.lstsq(outgoing, incoming, rcond=None)[0]
+
+
+@pytest.mark.parametrize('overlap', [.5, .8, .95])
+def test_truncated_mixed_polarizations_show_orientation_bias_despite_unitarity(overlap):
+    angle = np.arccos(overlap)
+    left, right = _polarizations(0)[:1], _polarizations(angle)[:1]
+    result = _scattering(left, right)
+    reversed_result = _reorder_ports(_scattering(right, left))
+    reflection = (1-overlap**2) / (1+overlap**2)
+    transmission = 2*overlap / (1+overlap**2)
+    assert_allclose(result, [[reflection, transmission], [transmission, -reflection]], atol=1e-14)
+    assert_allclose(reversed_result, [[-reflection, transmission], [transmission, reflection]], atol=1e-14)
+    assert_allclose(result.T @ result, np.eye(2), atol=1e-14)
+    assert_allclose(result, result.T, atol=1e-14)
+    assert_allclose(np.linalg.norm(result-reversed_result), 2*np.sqrt(2)*reflection, atol=1e-14)
+
+
+@pytest.mark.parametrize(('left_scales', 'right_scales'), [((1., 1.), (1., 1.)), ((1., 2.), (1., 3.))])
+def test_complete_maxwell_bases_recover_fresnel_and_port_relabeling(left_scales, right_scales):
+    left_angle, right_angle = np.pi/6, np.arccos(.8)
+    left, right = _polarizations(left_angle, left_scales), _polarizations(right_angle, right_scales)
+    # Establish that the synthetic fields are actual degenerate Maxwell modes,
+    # not arbitrary profiles or mixed vectors with different eigenvalues.
+    for modes, scales in ((left, left_scales), (right, right_scales)):
+        sx, sy = scales
+        epsilon, mu = np.array([1/sx**2, 1/sy**2, 1]), np.array([sy**2, sx**2, 1])
+        operator = waveguide_2d.operator_e(1, DXES, epsilon, mu)
+        for e, h in modes:
+            assert_allclose(operator @ e[:2], e[:2], atol=1e-14)
+            assert_allclose(waveguide_2d.exy2h(1, 1, DXES, epsilon, mu) @ e[:2], h, atol=1e-14)
+        gram = [[waveguide_2d.lorentz_overlap(a, b, DXES) for b in modes] for a in modes]
+        assert_allclose(gram, np.eye(2), atol=1e-14)
+    result = _scattering(left, right)
+    expected = _fresnel_reference(left_angle, right_angle, left_scales, right_scales)
+    assert_allclose(result, expected, atol=1e-14)
+    assert_allclose(_reorder_ports(_scattering(right, left)), expected, atol=1e-14)
+    assert_allclose(result.T @ result, np.eye(4), atol=1e-14)
+    assert_allclose(result, result.T, atol=1e-14)
+
+
+def test_missing_polarization_accounts_for_retained_submatrix_power_deficit():
+    overlap = .8
+    full = _scattering(_polarizations(0), _polarizations(np.arccos(overlap)))
+    retained = full[np.ix_([0, 2], [0, 2])]
+    assert_allclose(retained, [[0, overlap], [overlap, 0]], atol=1e-14)
+    assert_allclose(abs(full[3, 0])**2, 1-overlap**2, atol=1e-14)
+    assert_allclose(np.sum(abs(full[:, 0])**2), 1, atol=1e-14)
+    assert_allclose(np.sum(abs(retained[:, 0])**2), overlap**2, atol=1e-14)
+
+
+def test_posthoc_orientation_averaging_does_not_recover_physical_transmission():
+    overlap = .8
+    left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1]
+    averaged = .5 * (_scattering(left, right) + _reorder_ports(_scattering(right, left)))
+    assert_allclose(averaged[0, 0], 0, atol=1e-14)
+    assert_allclose(averaged[1, 1], 0, atol=1e-14)
+    assert abs(averaged[1, 0] - overlap) > .17
+
+
+def test_galerkin_testing_does_not_resolve_truncated_polarization_error():
+    overlap = .8
+    left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1]
+    result = _scattering(left, right, 'galerkin')
+    assert_allclose(result, [[0, overlap], [1/overlap, 0]], atol=1e-14)
+    assert not np.allclose(result, result.T)
+
+
+def test_trace_least_squares_recovers_uniform_reference_but_is_not_a_general_fix():
+    overlap = .8
+    left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1]
+    assert_allclose(_trace_least_squares(left, right), [[0, overlap], [overlap, 0]], atol=1e-14)
+    angle = np.pi/6
+    left, right = _polarizations(angle, (1., 2.)), _polarizations(angle, (1., 3.))
+    reference = _fresnel_reference(angle, angle, (1., 2.), (1., 3.))
+    assert_allclose(_trace_least_squares(left, right), reference, atol=1e-14)
+    truncated = _trace_least_squares(left[:1], right[:1])
+    # All traces come from the exact reciprocal, lossless Maxwell bases above.
+    # Euclidean residual minimization nevertheless violates both constraints
+    # after truncation. Port-relabeling covariance alone is not sufficient.
+    assert np.linalg.norm(truncated-truncated.T) > .3
+    assert np.linalg.svd(truncated, compute_uv=False).max() > 1.10
+    assert_allclose(_reorder_ports(_trace_least_squares(right[:1], left[:1])), truncated, atol=1e-14)
+
+
+@pytest.mark.parametrize('side', [0, 1])
+def test_driven_fdfd_confirms_uniform_polarization_reference(side):
+    omega, thickness = np.pi/8, 24
+    size = 2*thickness + 64
+    physical = [[np.ones(1), np.ones(1), np.ones(size)] for _ in range(2)]
+    epsilon = np.ones((3, 1, 1, size))
+    solve_widths = physical
+    for polarity in (-1, 1):
+        solve_widths = scpml.stretch_with_scpml(solve_widths, 2, polarity, omega, thickness=thickness)
+    overlap = .8
+    profiles = [np.array([1., 0.]), np.array([overlap, np.sqrt(1-overlap**2)])]
+    polarity = 1 if side == 0 else -1
+    source_index = thickness+4 if side == 0 else thickness+52
+    slices = (slice(None), slice(None), slice(source_index, source_index+1))
+    mode = waveguide_3d.mode_from_2d(profiles[side], omega, omega=omega, dxes=physical,
+                                    axis=2, polarity=polarity, slices=slices, epsilon=epsilon)
+    current = waveguide_3d.compute_source(mode['E'], mode['wavenumber'], omega, physical, 2, polarity, slices, epsilon)
+    field = unvec(solvers.generic(omega, solve_widths, vec(current), vec(epsilon), matrix_solver=spsolve), epsilon.shape[1:])
+    left_index, right_index = thickness+12, thickness+40
+    beta = 2*np.arcsin(omega/2)
+    phase = np.exp(-1j*beta)
+
+    def waves(index):
+        e0, e1 = field[:2, 0, 0, index], field[:2, 0, 0, index+1]
+        forward = (e1 - e0/phase) / (phase - 1/phase)
+        return forward, e0-forward
+
+    fl, bl = waves(left_index)
+    fr, br = waves(right_index)
+    incoming = np.dot(profiles[side], fl if side == 0 else br)
+    reflected = np.dot(profiles[side], bl if side == 0 else fr) / incoming
+    transmitted = np.dot(profiles[1-side], fr if side == 0 else bl) / incoming
+    transmitted *= np.exp(1j*beta*(right_index-left_index))
+    assert abs(reflected) < 5e-8
+    assert_allclose(transmitted, overlap, atol=1e-12)
+
+
+def test_longitudinal_refinement_cannot_repair_an_omitted_degenerate_polarization():
+    angle = np.arccos(.8)
+    reflection_magnitudes = []
+    for count in (3, 9, 33):
+        complete = [eme.ModeSection(z, _polarizations(angle*z), [1., 1.]) for z in np.linspace(0, 1, count)]
+        truncated = [eme.ModeSection(section.z, section.modes[:1], [1.]) for section in complete]
+        full = eme.get_taper_s(complete, DXES)
+        expected = _fresnel_reference(0, angle).astype(complex) * np.exp(-1j)
+        assert_allclose(full, expected, atol=1e-13)
+        small = eme.get_taper_s(truncated, DXES)
+        reflection_magnitudes.append(abs(small[0, 0]))
+        # There is no material taper: the basis is merely rotated within a
+        # degenerate subspace. Restricting the exact result gives |T|=.8, while
+        # the truncated chain approaches artificial adiabatic following, |T|=1.
+        assert abs(small[1, 0]) > .99
+    assert reflection_magnitudes[-1] < reflection_magnitudes[0] / 10
+    assert_allclose(abs(expected[2, 0]), .8, atol=1e-14)

meanas/test/test_examples_smoke.py

@@ -1,5 +1,6 @@
 from pathlib import Path
 import os
+import re
 import subprocess
 import sys
 
@@ -46,10 +47,14 @@ def test_eme_bend_example_smoke_runs(tmp_path: Path) -> None:
     assert 'straight effective indices:' in result.stdout
     assert 'cascaded bend through power' in result.stdout
 
+    straight_index_match = re.search(r'straight effective indices:\s*([0-9.]+)', result.stdout)
+    assert straight_index_match is not None, result.stdout
+    assert float(straight_index_match.group(1)) > 1.6
+
 
 def test_eme_taper_cmt_example_smoke_runs(tmp_path: Path) -> None:
     result = _run_example('eme_taper_cmt.py', tmp_path)
 
     assert result.returncode == 0, result.stdout + result.stderr
     assert 'sampled taper effective indices:' in result.stdout
-    assert 'taper CMT transmission' in result.stdout
+    assert 'taper transmission' in result.stdout

meanas/test/test_fdfd_algebra_helpers.py

@@ -51,10 +51,6 @@ def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
     return pm * (curl_term - OMEGA**2 * magnetic * masked_h)
 
 
-def _normalized_distance(u: numpy.ndarray, size: int, thickness: int) -> numpy.ndarray:
-    return ((thickness - u).clip(0) + (u - (size - thickness)).clip(0)) / thickness
-
-
 def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
     for mu in (None, MU):
         matrix_result = apply_fdfd_matrix(
@@ -155,14 +151,16 @@ def test_prepare_s_function_matches_closed_form_polynomial() -> None:
 def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None:
     s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
     dxes = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0, epsilon_effective=4.0, s_function=s_function)
-    correction = numpy.sqrt(4.0) * 2.0
-
     for axis, size, thickness in ((0, 6, 2), (2, 3, 1)):
-        grid = numpy.arange(size, dtype=float)
+        correction = thickness * numpy.sqrt(4.0) * 2.0
-        expected_a = 1 + 1j * s_function(_normalized_distance(grid, size, thickness)) / correction
+        for family, offset in enumerate((0.5, 0.0)):
-        expected_b = 1 + 1j * s_function(_normalized_distance(grid + 0.5, size, thickness)) / correction
+            positions = numpy.arange(size, dtype=float) + offset
-        assert_close(dxes[0][axis], expected_a)
+            expected = numpy.ones(size, dtype=complex)
-        assert_close(dxes[1][axis], expected_b)
+            low_depth = (thickness - positions[:thickness]) / thickness
+            high_depth = (positions[-thickness:] - (size - thickness)) / thickness
+            expected[:thickness] += 1j * s_function(low_depth) / correction
+            expected[-thickness:] += 1j * s_function(high_depth) / correction
+            assert_close(dxes[family][axis], expected)
 
     assert_close(dxes[0][1], 1.0)
     assert_close(dxes[1][1], 1.0)
@@ -182,7 +180,7 @@ def test_uniform_grid_scpml_default_s_function_matches_explicit_default() -> Non
 def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None:
     s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
     base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
-    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
+    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
 
     assert_close(stretched[0][0][2:], 1.0)
     assert_close(stretched[1][0][2:], 1.0)
@@ -195,7 +193,7 @@ def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None:
 def test_stretch_with_scpml_only_modifies_requested_back_edge() -> None:
     s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
     base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
-    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
+    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
 
     assert_close(stretched[0][0][:4], 1.0)
     assert_close(stretched[1][0][:4], 1.0)

meanas/test/test_numerical_feedback.py

@@ -0,0 +1,409 @@
+"""Independent regressions for the numerical/API feedback audit."""
+import dataclasses
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+from scipy import linalg, sparse
+
+from .. import eigensolvers
+from ..fdmath import vec
+from ..fdmath.functional import deriv_forward, deriv_back
+from ..fdtd import pml, base
+from ..fdfd import eme, scpml, waveguide_2d as wg, waveguide_3d
+
+
+def grid(shape):
+    return [[np.ones(n) for n in shape] for _ in range(2)]
+
+
+def mode(profile, admittance=1, scale=1):
+    profile = np.asarray(profile)
+    e = np.zeros((3, profile.size, 1), complex)
+    h = np.zeros_like(e)
+    e[0, :, 0] = np.sqrt(2) * scale * profile
+    h[1, :, 0] = np.sqrt(2) * scale * admittance * profile
+    return vec(e), vec(h)
+
+
+def test_disabled_cpml_matches_unequal_nonuniform_yee_updates():
+    rng = np.random.default_rng(6)
+    shape = (3, 4, 5, 6)
+    widths = [[np.linspace(.7 + .3 * family, 1.4 + .2 * family, n)
+               for n in shape[1:]] for family in range(2)]
+    e = rng.normal(size=shape) + 1j * rng.normal(size=shape)
+    h = rng.normal(size=shape) + 1j * rng.normal(size=shape)
+    eps, mu = 1 + rng.random(shape), 1 + rng.random(shape)
+    ec, hc = e.copy(), h.copy()
+    ue, uh = pml.updates_with_cpml([[None, None] for _ in range(3)], .1, widths, eps, dtype=complex)
+    ue(ec, hc, eps)
+    base.maxwell_e(.1, widths)(e, h, eps)
+    uh(ec, hc, mu)
+    base.maxwell_h(.1, widths)(e, h, mu)
+    assert_allclose(ec, e, atol=1e-14)
+    assert_allclose(hc, h, atol=1e-14)
+
+
+def test_cpml_physical_spacing_scales_attenuation_rate():
+    small = pml.cpml_params(0, 1, dt=.1, dx=.25)
+    unit = pml.cpml_params(0, 1, dt=.4)
+    for family in ('param_e', 'param_h'):
+        assert_allclose(small[family], unit[family])
+    for invalid in (0, -1, np.inf, np.nan):
+        with pytest.raises(ValueError, match='dx'):
+            pml.cpml_params(0, 1, .1, dx=invalid)
+
+
+def test_stretched_metric_satisfies_complex_reciprocity():
+    rng = np.random.default_rng(8)
+    shape = (4, 5)
+    widths = [[.7 + rng.random(n) - .2j * rng.random(n) for n in shape] for _ in range(2)]
+    operator = wg.operator_beta(.8, widths, vec(1 + rng.random((3, *shape))))
+    metric = wg.lorentz_metric(widths)
+    defect = operator.T @ metric - metric @ operator
+    assert sparse.linalg.norm(defect) / sparse.linalg.norm(metric @ operator) < 1e-12
+    physical = wg.lorentz_metric(widths, measure='physical')
+    assert sparse.linalg.norm(operator.T @ physical - physical @ operator) > 1e-3
+    e, h = mode([1, .3j, .2, .4])
+    # Direct one-dimensional complex-area quadrature, independent of the metric.
+    dx = [[np.full(4, 1 - .2j), np.ones(1)] for _ in range(2)]
+    assert_allclose(wg.inner_product(e, h, dx, measure='stretched'), (1 - .2j) * np.sum(np.array([1, .3j, .2, .4]) ** 2))
+    assert_allclose(wg.inner_product(e, h, dx, conj_h=True), np.sum(abs(np.array([1, .3j, .2, .4])) ** 2))
+
+
+@pytest.mark.parametrize('magnetic', [False, True])
+def test_reconstruction_can_normalize_complex_self_reaction(magnetic):
+    widths = [[np.full(2, 1 - .2j), np.ones(2)] for _ in range(2)]
+    transverse = np.concatenate((np.zeros(4), np.ones(4))) if magnetic else np.concatenate((np.ones(4), np.zeros(4)))
+    reconstruct = wg.normalized_fields_h if magnetic else wg.normalized_fields_e
+    e, h = reconstruct(transverse, 1, 1, widths, np.ones(12), normalization='lorentz')
+    assert_allclose(wg.inner_product(e, h, widths, measure='stretched'), 1, atol=1e-13)
+    with pytest.raises(ValueError, match='finite, nonzero'):
+        reconstruct(np.zeros(8), 1, 1, widths, np.ones(12))
+    with pytest.raises(ValueError, match='positive forward power'):
+        reconstruct(transverse, -1, 1, widths, np.ones(12))
+
+
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_scpml_face_copy_real_profile_and_final_stretch_cap(polarity):
+    widths = grid((10, 2))
+    result = scpml.stretch_with_scpml(
+        widths, 0, polarity, .2, thickness=3,
+        kappa_function=lambda depth: 1 + depth ** 2, max_imag_ratio=1,
+    )
+    selected = slice(0, 3) if polarity == -1 else slice(-3, None)
+    untouched = slice(3, None) if polarity == -1 else slice(None, -3)
+    assert result is not widths
+    for family in range(2):
+        assert_allclose(widths[family][0], 1)
+        assert_allclose(result[family][0][untouched], 1)
+        assert np.any(result[family][0][selected].real > 1)
+        assert np.any(result[family][0][selected].imag < 0)
+        assert np.all(abs(result[family][0].imag) <= result[family][0].real)
+        assert not np.shares_memory(result[family][1], widths[family][1])
+    mutated = scpml.stretch_with_scpml(widths, 0, polarity, .2, thickness=3, copy=False)
+    assert mutated is widths
+    assert np.any(widths[0][0].imag < 0)
+
+
+def test_uniform_scpml_caps_final_dimensionless_stretch():
+    widths = scpml.uniform_grid_scpml((10, 1, 1), (3, 0, 0), .2,
+                                     kappa_function=lambda depth: 1 + depth, max_imag_ratio=.7)
+    for family in widths:
+        assert np.all(abs(family[0].imag) <= .7 * family[0].real + 1e-14)
+        assert_allclose(family[1], 1)
+
+
+def test_indexed_eigensolve_accepts_seed_and_start_without_mutation():
+    operator = sparse.diags_array(np.linspace(-3, 5, 12))
+    first = eigensolvers.signed_eigensolve(operator, 2, rng=42)
+    second = eigensolvers.signed_eigensolve(operator, 2, rng=42)
+    assert_allclose(first[0], second[0], atol=1e-13)
+    assert_allclose(first[1], second[1], atol=1e-12)
+    start = np.linspace(1, 2, 12).astype(complex)
+    before = start.copy()
+    eigensolvers.signed_eigensolve(operator, 2, v0=start)
+    eigensolvers.power_iteration(operator, guess_vector=start)
+    assert_allclose(start, before)
+    with pytest.raises(ValueError, match='finite, nonzero'):
+        eigensolvers.signed_eigensolve(operator, 2, v0=np.zeros(12))
+
+
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_mode_from_2d_reuses_conversion_without_calling_eigensolver(monkeypatch, axis, polarity):
+    shape = (4, 5, 6)
+    widths, eps = grid(shape), np.ones((3, *shape))
+    slices = [slice(None)] * 3
+    slices[axis] = slice(1, 2)
+    points = int(np.prod([n for i, n in enumerate(shape) if i != axis]))
+    transverse = np.concatenate((np.ones(points), np.zeros(points)))
+    def forbidden(*args, **kwargs):
+        raise AssertionError('the conversion helper must not solve modes')
+    monkeypatch.setattr(wg, 'solve_mode', forbidden)
+    result = waveguide_3d.mode_from_2d(transverse, .5, omega=.5, dxes=widths, axis=axis,
+                                     polarity=polarity, slices=slices, epsilon=eps)
+    assert_allclose(result['wavenumber'], 2 * np.arcsin(.25))
+    assert np.isfinite(result['E']).all()
+    assert np.isfinite(result['H']).all()
+    monkeypatch.setattr(wg, 'solve_mode', lambda *args, **kwargs: (transverse, .5))
+    wrapped = waveguide_3d.solve_mode(0, .5, widths, axis, polarity, slices, eps)
+    assert_allclose(result['E'], wrapped['E'])
+    assert_allclose(result['H'], wrapped['H'])
+
+
+def test_mixed_projection_scalar_overlap_and_unclipped_galerkin():
+    a = .8
+    left, right = mode([1, 0]), mode([a, np.sqrt(1 - a*a)])
+    ss = eme.get_s([left], [1], [right], [1], grid((2, 1)))
+    transmission, reflection = 2*a/(1+a*a), (1-a*a)/(1+a*a)
+    assert_allclose(ss, [[reflection, transmission], [transmission, -reflection]], atol=1e-14)
+    assert_allclose(ss.T @ ss, np.eye(2), atol=1e-14)
+    t, r = eme.get_tr([left], [1], [right], [1], grid((2, 1)), formulation='galerkin')
+    assert_allclose(t, 1/a)
+    assert_allclose(r, 0, atol=1e-14)
+
+
+def test_asymmetric_transmission_uses_outgoing_rows_incoming_columns():
+    left, right = mode([1], scale=2), mode([1], scale=3)
+    ss = eme.get_s([left], [1], [right], [1], grid((1, 1)))
+    assert_allclose(ss, [[0, 1.5], [2/3, 0]], atol=1e-14)
+    # Continuity for the two physical incidence directions independently.
+    assert_allclose(2 * (1 + ss[0, 0]), 3 * ss[1, 0])
+    assert_allclose(3 * (1 + ss[1, 1]), 2 * ss[0, 1])
+    transfer = eme.get_abcd([left], [1], [right], [1], grid((1, 1))).toarray()
+    assert_allclose(transfer, np.eye(2) * (2/3), atol=1e-14)
+    assert_allclose(eme.abcd_to_s(transfer), ss, atol=1e-14)
+    assert_allclose(eme.abcd_to_s(sparse.csr_array(transfer)), ss, atol=1e-14)
+
+
+def test_full_nondiagonal_gram_projections_match_exact_field_coordinates():
+    left_matrix = np.array([[1, .2j], [.1, 1]], complex)
+    right_matrix = np.array([[1, .3], [.2, 1]], complex)
+    dl = np.array([[1, .2], [.3j, 1]], complex)
+    dr = np.array([[1, -.1j], [.4, 1]], complex)
+    left = [mode(column) for column in left_matrix.T]
+    right = [mode(column) for column in right_matrix.T]
+    ss = eme.get_s(left, [1, 1], right, [1, 1], grid((2, 1)),
+                   dual_ehLs=[mode(column) for column in dl.T],
+                   dual_ehRs=[mode(column) for column in dr.T])
+    zero = np.zeros((2, 2))
+    expected = np.block([[zero, np.linalg.solve(left_matrix, right_matrix)],
+                         [np.linalg.solve(right_matrix, left_matrix), zero]])
+    assert_allclose(ss, expected, atol=1e-13)
+
+
+def test_mixed_interface_accepts_unequal_port_counts():
+    left = [mode([1, 0])]
+    right = [mode([1, 0]), mode([0, 1])]
+    ss = eme.get_s(left, [1], right, [1, 2], grid((2, 1)))
+    assert_allclose(ss, [[0, 1, 0], [1, 0, 0], [0, 0, -1]], atol=1e-14)
+    with pytest.raises(ValueError, match='equal port mode counts'):
+        eme.get_abcd(left, [1], right, [1, 2], grid((2, 1)))
+
+
+def test_redheffer_matches_independent_internal_wave_elimination():
+    rng = np.random.default_rng(52)
+    left = .1 * (rng.normal(size=(3, 3)) + 1j * rng.normal(size=(3, 3)))
+    right = .1 * (rng.normal(size=(5, 5)) + 1j * rng.normal(size=(5, 5)))
+    internal = np.block([[np.eye(2), -left[1:, 1:]], [-right[:2, :2], np.eye(2)]])
+    incoming = np.block([[left[1:, :1], np.zeros((2, 3))], [np.zeros((2, 1)), right[:2, 2:]]])
+    waves = np.linalg.solve(internal, incoming)
+    outgoing = np.block([[left[:1, :1], np.zeros((1, 3))], [np.zeros((3, 1)), right[2:, 2:]]])
+    outgoing[:1] += left[:1, 1:] @ waves[2:]
+    outgoing[1:] += right[2:, :2] @ waves[:2]
+    assert_allclose(eme.cascade_s(left, right, internal_modes=2), outgoing, atol=1e-14)
+
+
+def test_taper_overlap_has_no_generator_gain():
+    left, right = mode([1, 0]), mode([.8, .6])
+    sections = [eme.ModeSection(0, [left], [2]), eme.ModeSection(1, [right], [2])]
+    ss = eme.get_taper_s(sections, grid((2, 1)))
+    assert_allclose(np.linalg.svd(ss, compute_uv=False), 1, atol=1e-14)
+    assert_allclose(abs(ss[1, 0]), 2*.8/(1+.8**2), atol=1e-14)
+    with pytest.raises(TypeError, match='rtol'):
+        eme.get_taper_s(sections, grid((2, 1)), rtol=1e-8)
+
+
+def test_long_leaky_taper_stays_finite_without_transfer_growth():
+    fields = mode([1])
+    sections = [eme.ModeSection(0, [fields], [1 - 1j]), eme.ModeSection(1000, [fields], [1 - 1j])]
+    ss = eme.get_taper_s(sections, grid((1, 1)))
+    assert_allclose(ss, 0, atol=1e-300)
+    with pytest.raises(np.linalg.LinAlgError, match='S to ABCD'):
+        eme.get_taper_abcd(sections, grid((1, 1)))
+
+
+def test_prepared_dense_beta_overrides_are_aligned_and_avoid_overlaps(monkeypatch):
+    fields = [mode([1, 0]), mode([0, 1])]
+    phases = np.exp(1j * np.array([.3, -.7]))
+    shifted = [(e * phase, h * phase) for (e, h), phase in zip(fields, phases, strict=True)]
+    beta = np.array([[1 - .1j, .2], [.1j, 2 - .1j]])
+    sections = [eme.ModeSection(0, fields, [1, 2], wavenumber_matrix=beta),
+                eme.ModeSection(2, shifted, [1, 2], wavenumber_matrix=beta)]
+    prepared = eme.prepare_taper(sections, grid((2, 1)))
+    overrides = np.array([beta * 1.1, beta * .9])
+    changed = [dataclasses.replace(section, z=z, wavenumber_matrix=b)
+               for section, z, b in zip(sections, [0, 3], overrides, strict=True)]
+    expected = eme.get_taper_s(changed, grid((2, 1)))
+    def forbidden(*args, **kwargs):
+        raise AssertionError('evaluation must not recompute modal overlaps')
+    monkeypatch.setattr(eme, 'inner_product', forbidden)
+    monkeypatch.setattr(eme, 'lorentz_overlap', forbidden)
+    actual = eme.evaluate_taper_s(prepared, z_coords=[0, 3], wavenumber_matrices=overrides)
+    assert_allclose(actual, expected, atol=1e-13)
+    fields[0][0][:] = 100
+    assert_allclose(eme.evaluate_taper_s(prepared, z_coords=[0, 3], wavenumber_matrices=overrides), expected, atol=1e-13)
+    with pytest.raises(ValueError, match='read-only'):
+        prepared.interfaces[0, 0, 0] = 5
+
+
+def test_explicit_backward_override_and_dense_uniform_propagation():
+    fields = [mode([1, 0]), mode([0, 1])]
+    backward = [(e, -h) for e, h in fields]
+    beta = np.array([[1 - .1j, .2], [.1j, 2 - .1j]])
+    minus = -beta + .03j * np.eye(2)
+    sections = [eme.ModeSection(z, fields, [1, 2], backward_modes=backward,
+                                backward_wavenumbers=[-1, -2], wavenumber_matrix=beta,
+                                backward_wavenumber_matrix=minus) for z in (0, 2)]
+    prepared = eme.prepare_taper(sections, grid((2, 1)))
+    with pytest.raises(ValueError, match='backward overrides'):
+        eme.evaluate_taper_s(prepared, wavenumber_matrices=np.array([beta, beta]))
+    ss = eme.evaluate_taper_s(prepared, wavenumber_matrices=np.array([beta, beta]),
+                              backward_wavenumber_matrices=np.array([minus, minus]))
+    expected = np.block([[np.zeros((2, 2)), linalg.expm(2j * minus)],
+                         [linalg.expm(-2j * beta), np.zeros((2, 2))]])
+    assert_allclose(ss, expected, atol=1e-13)
+
+
+@pytest.mark.parametrize('dtype', [np.int8, np.float32, np.float64, np.complex64, np.complex128])
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('forward', [False, True])
+@pytest.mark.parametrize('spaced', [False, True])
+@pytest.mark.parametrize('contiguous', [False, True])
+def test_slice_derivatives_match_periodic_reference_and_output_buffers(dtype, axis, forward, spaced, contiguous):
+    f = np.arange(6 * 8 * 4).reshape(6, 8, 4).astype(dtype)[::2, ::2, ::2]
+    if np.issubdtype(dtype, np.complexfloating):
+        f = f * (1 + .3j)
+    if contiguous:
+        f = np.ascontiguousarray(f)
+    widths = [np.linspace(.7, 1.3, n) for n in f.shape] if spaced else None
+    derivative = (deriv_forward if forward else deriv_back)(widths)[axis]
+    expected = np.roll(f, -1, axis=axis) - f if forward else f - np.roll(f, 1, axis=axis)
+    if spaced:
+        shape = [1, 1, 1]
+        shape[axis] = f.shape[axis]
+        expected = expected / widths[axis].reshape(shape)
+    assert_allclose(derivative(f), expected, rtol=1e-6, atol=1e-6)
+    output = np.empty_like(expected)
+    assert derivative(f, out=output) is output
+    assert_allclose(output, expected, rtol=1e-6, atol=1e-6)
+    with pytest.raises(ValueError, match='alias'):
+        derivative(f, out=f)
+    singleton = np.ones((1, 1, 1), dtype=dtype)
+    plain = (deriv_forward if forward else deriv_back)()[axis]
+    assert_allclose(plain(singleton), 0)
+
+
+def test_pml_cap_removes_above_light_line_artifact_in_homogeneous_fixture():
+    omega = .4
+    strong = grid((24, 1))
+    capped = grid((24, 1))
+    for polarity in (-1, 1):
+        strong = scpml.stretch_with_scpml(strong, 0, polarity, omega, thickness=5)
+        capped = scpml.stretch_with_scpml(capped, 0, polarity, omega, thickness=5, max_imag_ratio=1)
+    strong_values = np.linalg.eigvals(wg.operator_e(omega, strong, np.ones(72)).toarray())
+    capped_values = np.linalg.eigvals(wg.operator_e(omega, capped, np.ones(72)).toarray())
+    assert strong_values.real.max() > omega**2 + 1e-4
+    assert capped_values.real.max() <= omega**2 + 1e-10
+
+
+def test_sampled_taper_converges_to_continuous_maxwell_propagation():
+    from scipy.integrate import solve_ivp
+
+    omega, length = 1.3, 2.0
+    def maxwell(z, state):
+        index = 1 + z / length
+        operator = np.array([[0, -1j * omega], [-1j * omega * index**2, 0]])
+        return (operator @ state.reshape(2, 2)).ravel()
+    integrated = solve_ivp(maxwell, (0, length), np.eye(2, dtype=complex).ravel(), rtol=1e-11, atol=1e-13)
+    left_fields = np.array([[1, 1], [1, -1]], complex)
+    right_fields = np.array([[1, 1], [2, -2]], complex) / np.sqrt(2)
+    transfer = np.linalg.solve(right_fields, integrated.y[:, -1].reshape(2, 2) @ left_fields)
+    a, b, c, d = transfer.ravel()
+    reference = np.array([[-c/d, 1/d], [a-b*c/d, b/d]])
+    errors = []
+    for count in (9, 17, 33):
+        sections = []
+        for z in np.linspace(0, length, count):
+            index = 1 + z / length
+            fields = mode([1], admittance=index, scale=1/np.sqrt(index))
+            sections.append(eme.ModeSection(z, [fields], [omega * index]))
+        scattering = eme.get_taper_s(sections, grid((1, 1)))
+        assert_allclose(scattering.conj().T @ scattering, np.eye(2), atol=1e-12)
+        errors.append(np.linalg.norm(scattering - reference))
+    assert errors[1] < errors[0] / 3
+    assert errors[2] < errors[1] / 3
+    assert errors[2] < 5e-4
+
+
+def test_singular_conversions_and_interface_report_errors():
+    with pytest.raises(np.linalg.LinAlgError, match='S to ABCD'):
+        eme.s_to_abcd(np.eye(2))
+    with pytest.raises(np.linalg.LinAlgError, match='ABCD to S'):
+        eme.abcd_to_s(np.array([[1, 0], [0, 0]]))
+    fields = mode([0])
+    with pytest.raises(np.linalg.LinAlgError, match='interface continuity'):
+        eme.get_s([fields], [1], [fields], [1], grid((1, 1)))
+    mirror = np.eye(2)
+    with pytest.raises(np.linalg.LinAlgError, match='Redheffer internal feedback'):
+        eme.cascade_s(mirror, mirror, internal_modes=1)
+
+
+def test_uniform_custom_profiles_leave_physical_aperture_unchanged():
+    widths = scpml.uniform_grid_scpml((10, 1, 1), (2, 0, 0), 1,
+                                     kappa_function=lambda depth: np.full_like(depth, 2))
+    for family in widths:
+        assert_allclose(family[0][3:8], 1)
+
+
+@pytest.mark.parametrize('forward', [False, True])
+def test_derivatives_do_not_subtract_across_independent_rows(forward):
+    fields = np.array([[[3e38, 3e38], [-3e38, -3e38]]], dtype=np.float32)
+    derivative = (deriv_forward if forward else deriv_back)()[2]
+    with np.errstate(over='raise', invalid='raise'):
+        assert_allclose(derivative(fields), 0)
+
+
+@pytest.mark.parametrize(('shape', 'thicknesses'), [
+    ((6, 8, 4), (0, 0, 0)),
+    ((6, 8, 4), (1, 1, 1)),
+    ((6, 8, 4), (2, 3, 0)),
+    ((6, 8, 4), (3, 4, 2)),
+])
+@pytest.mark.parametrize('kappa', [None, 'polynomial', 'constant'])
+@pytest.mark.parametrize('cap', [None, .7, 1.0])
+def test_scpml_builders_agree_including_real_profiles_and_touching_faces(shape, thicknesses, kappa, cap):
+    profile = scpml.prepare_s_function(ln_R=-9, m=2.5)
+    real_profile = {
+        None: None,
+        'polynomial': lambda depth: 1 + 2 * depth**2,
+        'constant': lambda depth: np.full_like(depth, 2),
+    }[kappa]
+    options = dict(omega=.8, epsilon_effective=2.25, s_function=profile,
+                   kappa_function=real_profile, max_imag_ratio=cap)
+    uniform = scpml.uniform_grid_scpml(shape, thicknesses, **options)
+    explicit = grid(shape)
+    for axis, thickness in enumerate(thicknesses):
+        for polarity in (-1, 1):
+            explicit = scpml.stretch_with_scpml(explicit, axis, polarity, thickness=thickness, **options)
+    for family in range(2):
+        for axis in range(3):
+            assert_allclose(uniform[family][axis], explicit[family][axis], rtol=1e-14, atol=1e-14)
+
+
+@pytest.mark.parametrize('attenuation', [0.0, -0.0, .1, np.inf, -np.inf, np.nan])
+@pytest.mark.parametrize('alpha', [0, .2])
+def test_cpml_rejects_zero_or_invalid_attenuation_before_coefficient_evaluation(attenuation, alpha):
+    with np.errstate(divide='raise', invalid='raise'), pytest.raises(ValueError, match='ln_R_per_layer must be finite and negative'):
+        pml.cpml_params(0, 1, .1, ln_R_per_layer=attenuation, cfs_alpha=alpha)

meanas/test/test_waveguide_2d_numerics.py

@@ -129,6 +129,7 @@ def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu
         omega=OMEGA,
         dxes=DXES_2D,
         epsilon=epsilon,
+        mu=mu,
         )
     _e_ref, h_ref = waveguide_2d.normalized_fields_e(
         e_xy,
@@ -167,6 +168,7 @@ def test_waveguide_2d_helper_operators_with_mu_return_finite_outputs() -> None:
         omega=OMEGA,
         dxes=DXES_2D,
         epsilon=epsilon,
+        mu=mu,
         )
     _e_ref, h_ref = waveguide_2d.normalized_fields_e(
         e_xy,
@@ -202,6 +204,7 @@ def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None:
         omega=OMEGA,
         dxes=DXES_2D,
         epsilon=epsilon,
+        mu=mu,
         )
     e_field, h_field = waveguide_2d.normalized_fields_e(
         e_xy,
@@ -217,8 +220,8 @@ def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None:
 
     assert numpy.isfinite(h_error)
     assert numpy.isfinite(e_error)
-    assert 0 < h_error < 0.1
+    assert 0 < h_error < 1e-6
-    assert 0 < e_error < 2.0
+    assert 0 < e_error < 1e-6
 
 
 def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None:
@@ -230,6 +233,7 @@ def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None:
         omega=OMEGA,
         dxes=DXES_2D,
         epsilon=epsilon,
+        mu=mu,
         )
     e_field, h_field = waveguide_2d.normalized_fields_e(
         e_xy,
@@ -282,3 +286,35 @@ def test_waveguide_2d_inner_product_trapezoid_branch_on_nonuniform_grid() -> Non
     assert abs(overlap_rect.imag) < 1e-10
     assert abs(overlap_trap.imag) < 1e-10
     assert abs(overlap_rect - overlap_trap) > 0.1
+
+
+def test_waveguide_2d_inner_product_uses_staggered_nonuniform_areas() -> None:
+    shape = (3, 4)
+    dxes = [
+        [numpy.asarray([0.4, 0.7, 1.1]), numpy.asarray([0.3, 0.5, 0.8, 1.2])],
+        [numpy.asarray([0.6, 0.9, 1.3]), numpy.asarray([0.2, 0.45, 0.75, 1.0])],
+        ]
+    samples = numpy.arange(numpy.prod(shape), dtype=float).reshape(shape)
+    electric = vec(numpy.stack((1.0 + samples, 0.25 - 0.3 * samples, numpy.zeros(shape))))
+    magnetic = vec(numpy.stack((-0.4 + 0.2 * samples, 2.0 - 0.1 * samples, numpy.zeros(shape))))
+    ex, ey, _ez = numpy.split(electric, 3)
+    hx, hy, _hz = numpy.split(magnetic, 3)
+    ex = ex.reshape(shape)
+    ey = ey.reshape(shape)
+    hx = hx.reshape(shape)
+    hy = hy.reshape(shape)
+
+    expected_rectangular = 0.5 * (
+        numpy.sum(ex * hy * dxes[0][0][:, None] * dxes[1][1][None, :])
+        - numpy.sum(ey * hx * dxes[1][0][:, None] * dxes[0][1][None, :])
+        )
+    expected_trapezoidal = 0.5 * (
+        numpy.trapezoid(numpy.trapezoid(ex * hy, numpy.cumsum(dxes[1][1])), numpy.cumsum(dxes[0][0]))
+        - numpy.trapezoid(numpy.trapezoid(ey * hx, numpy.cumsum(dxes[0][1])), numpy.cumsum(dxes[1][0]))
+        )
+
+    assert_allclose(waveguide_2d.inner_product(electric, magnetic, dxes), expected_rectangular)
+    assert_allclose(
+        waveguide_2d.inner_product(electric, magnetic, dxes, trapezoid=True),
+        expected_trapezoidal,
+        )

meanas/test/test_waveguide_mode_helpers.py

@@ -6,6 +6,7 @@ import pytest
 import warnings
 
 from ..fdmath import vec, unvec
+from ..fdmath.operators import deriv_forward, deriv_back
 from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl
 
 
@@ -16,15 +17,21 @@ def build_waveguide_3d_mode(
         *,
         slice_start: int,
         polarity: int,
+        axis: int = 0,
+        shape: tuple[int, int, int] = (5, 5, 1),
+        dxes: list[list[numpy.ndarray]] | None = None,
         ) -> tuple[numpy.ndarray, list[list[numpy.ndarray]], tuple[slice, slice, slice], waveguide_3d.Waveguide3DMode]:
-    epsilon = numpy.ones((3, 5, 5, 1), dtype=float)
+    epsilon = numpy.ones((3, *shape), dtype=float)
-    dxes = [[numpy.ones(5), numpy.ones(5), numpy.ones(1)] for _ in range(2)]
+    if dxes is None:
-    slices = (slice(slice_start, slice_start + 1), slice(None), slice(None))
+        dxes = [[numpy.ones(size) for size in shape] for _ in range(2)]
+    slices_l = [slice(None), slice(None), slice(None)]
+    slices_l[axis] = slice(slice_start, slice_start + 1)
+    slices = tuple(slices_l)
     result = waveguide_3d.solve_mode(
         0,
         omega=OMEGA,
         dxes=dxes,
-        axis=0,
+        axis=axis,
         polarity=polarity,
         slices=slices,
         epsilon=epsilon,
@@ -55,10 +62,37 @@ def build_waveguide_cyl_mu_profile() -> numpy.ndarray:
     return numpy.linspace(1.5, 2.2, 3 * 5 * 5)
 
 
-def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
+def relative_residual(lhs: numpy.ndarray, rhs: numpy.ndarray) -> float:
-    epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1)
+    return float(norm(lhs - rhs) / (norm(lhs) + norm(rhs)))
-    axis = 0
+
+
+@pytest.mark.parametrize(
+    ('axis', 'shape'),
+    [
+        (0, (5, 5, 1)),
+        (1, (5, 5, 1)),
+        (2, (5, 5, 5)),
+    ],
+    )
+def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent(
+        axis: int,
+        shape: tuple[int, int, int],
+        ) -> None:
     polarity = 1
+    dxes = [
+        [
+            numpy.linspace(0.9 + 0.1 * eh + 0.05 * aa, 1.2 + 0.1 * eh + 0.05 * aa, size)
+            for aa, size in enumerate(shape)
+        ]
+        for eh in range(2)
+    ]
+    epsilon, dxes, slices, result = build_waveguide_3d_mode(
+        slice_start=2,
+        polarity=polarity,
+        axis=axis,
+        shape=shape,
+        dxes=dxes,
+        )
     expanded = waveguide_3d.expand_e(
         E=result['E'],
         wavenumber=result['wavenumber'],
@@ -68,7 +102,7 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
         slices=slices,
         )
 
-    dx_prop = 0.5 * numpy.array([dx[2][slices[2]] for dx in dxes]).sum()
+    dx_prop = 0.5 * numpy.array([dx[axis][slices[axis]] for dx in dxes]).sum()
     expected_wavenumber = 2 / dx_prop * numpy.arcsin(result['wavenumber_2d'] * dx_prop / 2)
     solved_slice = (slice(None), *slices)
 
@@ -79,15 +113,18 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
     assert abs(result['wavenumber'] - expected_wavenumber) < 1e-12
     assert numpy.allclose(expanded[solved_slice], result['E'][solved_slice])
 
-    component, _x, y_index, z_index = numpy.unravel_index(
+    component, *indices = numpy.unravel_index(
         numpy.abs(result['E']).argmax(),
         result['E'].shape,
         )
-    values = expanded[component, :, y_index, z_index]
+    line_selection: list[int | slice] = list(indices)
+    line_selection[axis] = slice(None)
+    values = expanded[component][tuple(line_selection)]
     ratios = values[1:] / values[:-1]
-    expected_ratio = numpy.exp(-1j * result['wavenumber'])
+    phase_positions = numpy.real(dxes[0][axis]).cumsum()
+    expected_ratios = numpy.exp(-1j * result['wavenumber'] * numpy.diff(phase_positions))
 
-    numpy.testing.assert_allclose(ratios, expected_ratio, rtol=1e-6, atol=1e-9)
+    numpy.testing.assert_allclose(ratios, expected_ratios, rtol=1e-6, atol=1e-9)
 
 
 @pytest.mark.parametrize(
@@ -117,6 +154,17 @@ def test_waveguide_3d_compute_overlap_e_uses_adjacent_window(
     assert nonzero[:, 1].min() == expected_range[0]
     assert nonzero[:, 1].max() == expected_range[1]
 
+    expanded = waveguide_3d.expand_e(
+        E=result['E'],
+        wavenumber=result['wavenumber'],
+        dxes=dxes,
+        axis=0,
+        polarity=polarity,
+        slices=slices,
+        )
+    numpy.testing.assert_allclose(numpy.sum(overlap * expanded.conj()), 1.0, atol=1e-12)
+    numpy.testing.assert_allclose(vec(expanded) @ vec(overlap).conj(), 1.0, atol=1e-12)
+
 
 @pytest.mark.parametrize(
     ('polarity', 'slice_start', 'expected_index'),
@@ -301,6 +349,61 @@ def test_waveguide_cyl_reconstruction_matches_2d_limit(use_mu: bool) -> None:
     numpy.testing.assert_allclose(cyl_h, straight_h, rtol=1e-8, atol=1e-8)
 
 
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_waveguide_cyl_reconstruction_is_finite_radius_consistent(use_mu: bool) -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True, asymmetric=True)
+    mu = build_waveguide_cyl_mu_profile() if use_mu else None
+    mu_full = numpy.ones_like(epsilon) if mu is None else mu
+
+    e_xy, angular_wavenumber = waveguide_cyl.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=rmin,
+        mu=mu,
+        )
+    e_field = waveguide_cyl.exy2e(
+        angular_wavenumber=angular_wavenumber,
+        dxes=dxes,
+        rmin=rmin,
+        epsilon=epsilon,
+        ) @ e_xy
+    h_field = waveguide_cyl.exy2h(
+        angular_wavenumber=angular_wavenumber,
+        omega=OMEGA,
+        dxes=dxes,
+        rmin=rmin,
+        epsilon=epsilon,
+        mu=mu,
+        ) @ e_xy
+
+    Dfx, Dfy = deriv_forward(dxes[0])
+    Dbx, Dby = deriv_back(dxes[1])
+    Ta, Tb = waveguide_cyl.dxes2T(dxes, rmin)
+    ta = Ta.diagonal()
+    tb = Tb.diagonal()
+    beta = angular_wavenumber / rmin
+
+    er, ey, ez = numpy.split(e_field, 3)
+    hr, hy, hz = numpy.split(h_field, 3)
+    eps_r, eps_y, eps_z = numpy.split(epsilon, 3)
+    mu_r, mu_y, mu_z = numpy.split(mu_full, 3)
+
+    ez_rhs = ta * (
+        (Dbx @ (tb * eps_r * er))
+        + (Dby @ (ta * eps_y * ey))
+        ) / eps_z
+    hr_rhs = 1j * beta * ey / ta + Dfy @ ez
+    hy_rhs = -1j * beta * er / tb - (Dfx @ (ta * ez)) / tb
+    hz_rhs = -Dfy @ er + Dfx @ ey
+
+    assert relative_residual(1j * beta * ez, ez_rhs) < 1e-7
+    assert relative_residual(-1j * OMEGA * mu_r * hr, hr_rhs) < 1e-7
+    assert relative_residual(-1j * OMEGA * mu_y * hy, hy_rhs) < 1e-7
+    assert relative_residual(-1j * OMEGA * mu_z * hz, hz_rhs) < 1e-7
+
+
 def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None:
     dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True)
     nx = dxes[0][0].size

pyproject.toml

@@ -33,7 +33,7 @@ classifiers = [
     "License :: OSI Approved :: GNU Affero General Public License v3",
     "Topic :: Scientific/Engineering :: Physics",
     ]
-requires-python = ">=3.11"
+requires-python = ">=3.12"
 include = [
     "LICENSE.md"
     ]

uv.lock

@@ -1,5 +1,5 @@
 version = 1
-requires-python = ">=3.11"
+requires-python = ">=3.12"
 resolution-markers = [
     "python_full_version >= '3.14' and sys_platform == 'win32'",
     "python_full_version >= '3.14' and sys_platform == 'emscripten'",
@@ -60,22 +60,6 @@ version = "3.4.7"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/e7/a1/67fe25fac3c7642725500a3f6cfe5821ad557c3abb11c9d20d12c7008d3e/charset_normalizer-3.4.7.tar.gz", hash = "sha256:ae89db9e5f98a11a4bf50407d4363e7b09b31e55bc117b4f7d80aab97ba009e5", size = 144271 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/c2/d7/b5b7020a0565c2e9fa8c09f4b5fa6232feb326b8c20081ccded47ea368fd/charset_normalizer-3.4.7-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:7641bb8895e77f921102f72833904dcd9901df5d6d72a2ab8f31d04b7e51e4e7", size = 309705 },
-    { url = "https://files.pythonhosted.org/packages/5a/53/58c29116c340e5456724ecd2fff4196d236b98f3da97b404bc5e51ac3493/charset_normalizer-3.4.7-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:202389074300232baeb53ae2569a60901f7efadd4245cf3a3bf0617d60b439d7", size = 206419 },
-    { url = "https://files.pythonhosted.org/packages/b2/02/e8146dc6591a37a00e5144c63f29fb7c97a734ea8a111190783c0e60ab63/charset_normalizer-3.4.7-cp311-cp311-manylinux2014_ppc64le.manylinux_2_17_ppc64le.manylinux_2_28_ppc64le.whl", hash = "sha256:30b8d1d8c52a48c2c5690e152c169b673487a2a58de1ec7393196753063fcd5e", size = 227901 },
-    { url = "https://files.pythonhosted.org/packages/fb/73/77486c4cd58f1267bf17db420e930c9afa1b3be3fe8c8b8ebbebc9624359/charset_normalizer-3.4.7-cp311-cp311-manylinux2014_s390x.manylinux_2_17_s390x.manylinux_2_28_s390x.whl", hash = "sha256:532bc9bf33a68613fd7d65e4b1c71a6a38d7d42604ecf239c77392e9b4e8998c", size = 222742 },
-    { url = "https://files.pythonhosted.org/packages/a1/fa/f74eb381a7d94ded44739e9d94de18dc5edc9c17fb8c11f0a6890696c0a9/charset_normalizer-3.4.7-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:2fe249cb4651fd12605b7288b24751d8bfd46d35f12a20b1ba33dea122e690df", size = 214061 },
-    { url = "https://files.pythonhosted.org/packages/dc/92/42bd3cefcf7687253fb86694b45f37b733c97f59af3724f356fa92b8c344/charset_normalizer-3.4.7-cp311-cp311-manylinux_2_31_armv7l.whl", hash = "sha256:65bcd23054beab4d166035cabbc868a09c1a49d1efe458fe8e4361215df40265", size = 199239 },
-    { url = "https://files.pythonhosted.org/packages/4c/3d/069e7184e2aa3b3cddc700e3dd267413dc259854adc3380421c805c6a17d/charset_normalizer-3.4.7-cp311-cp311-manylinux_2_31_riscv64.manylinux_2_39_riscv64.whl", hash = "sha256:08e721811161356f97b4059a9ba7bafb23ea5ee2255402c42881c214e173c6b4", size = 210173 },
-    { url = "https://files.pythonhosted.org/packages/62/51/9d56feb5f2e7074c46f93e0ebdbe61f0848ee246e2f0d89f8e20b89ebb8f/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:e060d01aec0a910bdccb8be71faf34e7799ce36950f8294c8bf612cba65a2c9e", size = 209841 },
-    { url = "https://files.pythonhosted.org/packages/d2/59/893d8f99cc4c837dda1fe2f1139079703deb9f321aabcb032355de13b6c7/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_armv7l.whl", hash = "sha256:38c0109396c4cfc574d502df99742a45c72c08eff0a36158b6f04000043dbf38", size = 200304 },
-    { url = "https://files.pythonhosted.org/packages/7d/1d/ee6f3be3464247578d1ed5c46de545ccc3d3ff933695395c402c21fa6b77/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_ppc64le.whl", hash = "sha256:1c2a768fdd44ee4a9339a9b0b130049139b8ce3c01d2ce09f67f5a68048d477c", size = 229455 },
-    { url = "https://files.pythonhosted.org/packages/54/bb/8fb0a946296ea96a488928bdce8ef99023998c48e4713af533e9bb98ef07/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_riscv64.whl", hash = "sha256:1a87ca9d5df6fe460483d9a5bbf2b18f620cbed41b432e2bddb686228282d10b", size = 210036 },
-    { url = "https://files.pythonhosted.org/packages/9a/bc/015b2387f913749f82afd4fcba07846d05b6d784dd16123cb66860e0237d/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_s390x.whl", hash = "sha256:d635aab80466bc95771bb78d5370e74d36d1fe31467b6b29b8b57b2a3cd7d22c", size = 224739 },
-    { url = "https://files.pythonhosted.org/packages/17/ab/63133691f56baae417493cba6b7c641571a2130eb7bceba6773367ab9ec5/charset_normalizer-3.4.7-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:ae196f021b5e7c78e918242d217db021ed2a6ace2bc6ae94c0fc596221c7f58d", size = 216277 },
-    { url = "https://files.pythonhosted.org/packages/06/6d/3be70e827977f20db77c12a97e6a9f973631a45b8d186c084527e53e77a4/charset_normalizer-3.4.7-cp311-cp311-win32.whl", hash = "sha256:adb2597b428735679446b46c8badf467b4ca5f5056aae4d51a19f9570301b1ad", size = 147819 },
-    { url = "https://files.pythonhosted.org/packages/20/d9/5f67790f06b735d7c7637171bbfd89882ad67201891b7275e51116ed8207/charset_normalizer-3.4.7-cp311-cp311-win_amd64.whl", hash = "sha256:8e385e4267ab76874ae30db04c627faaaf0b509e1ccc11a95b3fc3e83f855c00", size = 159281 },
-    { url = "https://files.pythonhosted.org/packages/ca/83/6413f36c5a34afead88ce6f66684d943d91f233d76dd083798f9602b75ae/charset_normalizer-3.4.7-cp311-cp311-win_arm64.whl", hash = "sha256:d4a48e5b3c2a489fae013b7589308a40146ee081f6f509e047e0e096084ceca1", size = 147843 },
     { url = "https://files.pythonhosted.org/packages/0c/eb/4fc8d0a7110eb5fc9cc161723a34a8a6c200ce3b4fbf681bc86feee22308/charset_normalizer-3.4.7-cp312-cp312-macosx_10_13_universal2.whl", hash = "sha256:eca9705049ad3c7345d574e3510665cb2cf844c2f2dcfe675332677f081cbd46", size = 311328 },
     { url = "https://files.pythonhosted.org/packages/f8/e3/0fadc706008ac9d7b9b5be6dc767c05f9d3e5df51744ce4cc9605de7b9f4/charset_normalizer-3.4.7-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:6178f72c5508bfc5fd446a5905e698c6212932f25bcdd4b47a757a50605a90e2", size = 208061 },
     { url = "https://files.pythonhosted.org/packages/42/f0/3dd1045c47f4a4604df85ec18ad093912ae1344ac706993aff91d38773a2/charset_normalizer-3.4.7-cp312-cp312-manylinux2014_ppc64le.manylinux_2_17_ppc64le.manylinux_2_28_ppc64le.whl", hash = "sha256:e1421b502d83040e6d7fb2fb18dff63957f720da3d77b2fbd3187ceb63755d7b", size = 229031 },
@@ -173,17 +157,6 @@ dependencies = [
 ]
 sdist = { url = "https://files.pythonhosted.org/packages/58/01/1253e6698a07380cd31a736d248a3f2a50a7c88779a1813da27503cadc2a/contourpy-1.3.3.tar.gz", hash = "sha256:083e12155b210502d0bca491432bb04d56dc3432f95a979b429f2848c3dbe880", size = 13466174 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/91/2e/c4390a31919d8a78b90e8ecf87cd4b4c4f05a5b48d05ec17db8e5404c6f4/contourpy-1.3.3-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:709a48ef9a690e1343202916450bc48b9e51c049b089c7f79a267b46cffcdaa1", size = 288773 },
-    { url = "https://files.pythonhosted.org/packages/0d/44/c4b0b6095fef4dc9c420e041799591e3b63e9619e3044f7f4f6c21c0ab24/contourpy-1.3.3-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:23416f38bfd74d5d28ab8429cc4d63fa67d5068bd711a85edb1c3fb0c3e2f381", size = 270149 },
-    { url = "https://files.pythonhosted.org/packages/30/2e/dd4ced42fefac8470661d7cb7e264808425e6c5d56d175291e93890cce09/contourpy-1.3.3-cp311-cp311-manylinux_2_26_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:929ddf8c4c7f348e4c0a5a3a714b5c8542ffaa8c22954862a46ca1813b667ee7", size = 329222 },
-    { url = "https://files.pythonhosted.org/packages/f2/74/cc6ec2548e3d276c71389ea4802a774b7aa3558223b7bade3f25787fafc2/contourpy-1.3.3-cp311-cp311-manylinux_2_26_ppc64le.manylinux_2_28_ppc64le.whl", hash = "sha256:9e999574eddae35f1312c2b4b717b7885d4edd6cb46700e04f7f02db454e67c1", size = 377234 },
-    { url = "https://files.pythonhosted.org/packages/03/b3/64ef723029f917410f75c09da54254c5f9ea90ef89b143ccadb09df14c15/contourpy-1.3.3-cp311-cp311-manylinux_2_26_s390x.manylinux_2_28_s390x.whl", hash = "sha256:0bf67e0e3f482cb69779dd3061b534eb35ac9b17f163d851e2a547d56dba0a3a", size = 380555 },
-    { url = "https://files.pythonhosted.org/packages/5f/4b/6157f24ca425b89fe2eb7e7be642375711ab671135be21e6faa100f7448c/contourpy-1.3.3-cp311-cp311-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:51e79c1f7470158e838808d4a996fa9bac72c498e93d8ebe5119bc1e6becb0db", size = 355238 },
-    { url = "https://files.pythonhosted.org/packages/98/56/f914f0dd678480708a04cfd2206e7c382533249bc5001eb9f58aa693e200/contourpy-1.3.3-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:598c3aaece21c503615fd59c92a3598b428b2f01bfb4b8ca9c4edeecc2438620", size = 1326218 },
-    { url = "https://files.pythonhosted.org/packages/fb/d7/4a972334a0c971acd5172389671113ae82aa7527073980c38d5868ff1161/contourpy-1.3.3-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:322ab1c99b008dad206d406bb61d014cf0174df491ae9d9d0fac6a6fda4f977f", size = 1392867 },
-    { url = "https://files.pythonhosted.org/packages/75/3e/f2cc6cd56dc8cff46b1a56232eabc6feea52720083ea71ab15523daab796/contourpy-1.3.3-cp311-cp311-win32.whl", hash = "sha256:fd907ae12cd483cd83e414b12941c632a969171bf90fc937d0c9f268a31cafff", size = 183677 },
-    { url = "https://files.pythonhosted.org/packages/98/4b/9bd370b004b5c9d8045c6c33cf65bae018b27aca550a3f657cdc99acdbd8/contourpy-1.3.3-cp311-cp311-win_amd64.whl", hash = "sha256:3519428f6be58431c56581f1694ba8e50626f2dd550af225f82fb5f5814d2a42", size = 225234 },
-    { url = "https://files.pythonhosted.org/packages/d9/b6/71771e02c2e004450c12b1120a5f488cad2e4d5b590b1af8bad060360fe4/contourpy-1.3.3-cp311-cp311-win_arm64.whl", hash = "sha256:15ff10bfada4bf92ec8b31c62bf7c1834c244019b4a33095a68000d7075df470", size = 193123 },
     { url = "https://files.pythonhosted.org/packages/be/45/adfee365d9ea3d853550b2e735f9d66366701c65db7855cd07621732ccfc/contourpy-1.3.3-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:b08a32ea2f8e42cf1d4be3169a98dd4be32bafe4f22b6c4cb4ba810fa9e5d2cb", size = 293419 },
     { url = "https://files.pythonhosted.org/packages/53/3e/405b59cfa13021a56bba395a6b3aca8cec012b45bf177b0eaf7a202cde2c/contourpy-1.3.3-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:556dba8fb6f5d8742f2923fe9457dbdd51e1049c4a43fd3986a0b14a1d815fc6", size = 273979 },
     { url = "https://files.pythonhosted.org/packages/d4/1c/a12359b9b2ca3a845e8f7f9ac08bdf776114eb931392fcad91743e2ea17b/contourpy-1.3.3-cp312-cp312-manylinux_2_26_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:92d9abc807cf7d0e047b95ca5d957cf4792fcd04e920ca70d48add15c1a90ea7", size = 332653 },
@@ -239,11 +212,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/93/8a/68a4ec5c55a2971213d29a9374913f7e9f18581945a7a31d1a39b5d2dfe5/contourpy-1.3.3-cp314-cp314t-win32.whl", hash = "sha256:e74a9a0f5e3fff48fb5a7f2fd2b9b70a3fe014a67522f79b7cca4c0c7e43c9ae", size = 202428 },
     { url = "https://files.pythonhosted.org/packages/fa/96/fd9f641ffedc4fa3ace923af73b9d07e869496c9cc7a459103e6e978992f/contourpy-1.3.3-cp314-cp314t-win_amd64.whl", hash = "sha256:13b68d6a62db8eafaebb8039218921399baf6e47bf85006fd8529f2a08ef33fc", size = 250331 },
     { url = "https://files.pythonhosted.org/packages/ae/8c/469afb6465b853afff216f9528ffda78a915ff880ed58813ba4faf4ba0b6/contourpy-1.3.3-cp314-cp314t-win_arm64.whl", hash = "sha256:b7448cb5a725bb1e35ce88771b86fba35ef418952474492cf7c764059933ff8b", size = 203831 },
-    { url = "https://files.pythonhosted.org/packages/a5/29/8dcfe16f0107943fa92388c23f6e05cff0ba58058c4c95b00280d4c75a14/contourpy-1.3.3-pp311-pypy311_pp73-macosx_10_15_x86_64.whl", hash = "sha256:cd5dfcaeb10f7b7f9dc8941717c6c2ade08f587be2226222c12b25f0483ed497", size = 278809 },
-    { url = "https://files.pythonhosted.org/packages/85/a9/8b37ef4f7dafeb335daee3c8254645ef5725be4d9c6aa70b50ec46ef2f7e/contourpy-1.3.3-pp311-pypy311_pp73-macosx_11_0_arm64.whl", hash = "sha256:0c1fc238306b35f246d61a1d416a627348b5cf0648648a031e14bb8705fcdfe8", size = 261593 },
-    { url = "https://files.pythonhosted.org/packages/0a/59/ebfb8c677c75605cc27f7122c90313fd2f375ff3c8d19a1694bda74aaa63/contourpy-1.3.3-pp311-pypy311_pp73-manylinux_2_26_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:70f9aad7de812d6541d29d2bbf8feb22ff7e1c299523db288004e3157ff4674e", size = 302202 },
-    { url = "https://files.pythonhosted.org/packages/3c/37/21972a15834d90bfbfb009b9d004779bd5a07a0ec0234e5ba8f64d5736f4/contourpy-1.3.3-pp311-pypy311_pp73-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:5ed3657edf08512fc3fe81b510e35c2012fbd3081d2e26160f27ca28affec989", size = 329207 },
-    { url = "https://files.pythonhosted.org/packages/0c/58/bd257695f39d05594ca4ad60df5bcb7e32247f9951fd09a9b8edb82d1daa/contourpy-1.3.3-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:3d1a3799d62d45c18bafd41c5fa05120b96a28079f2393af559b843d1a966a77", size = 225315 },
 ]
 
 [[package]]
@@ -252,21 +220,6 @@ version = "7.13.5"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/9d/e0/70553e3000e345daff267cec284ce4cbf3fc141b6da229ac52775b5428f1/coverage-7.13.5.tar.gz", hash = "sha256:c81f6515c4c40141f83f502b07bbfa5c240ba25bbe73da7b33f1e5b6120ff179", size = 915967 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/4b/37/d24c8f8220ff07b839b2c043ea4903a33b0f455abe673ae3c03bbdb7f212/coverage-7.13.5-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:66a80c616f80181f4d643b0f9e709d97bcea413ecd9631e1dedc7401c8e6695d", size = 219381 },
-    { url = "https://files.pythonhosted.org/packages/35/8b/cd129b0ca4afe886a6ce9d183c44d8301acbd4ef248622e7c49a23145605/coverage-7.13.5-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:145ede53ccbafb297c1c9287f788d1bc3efd6c900da23bf6931b09eafc931587", size = 219880 },
-    { url = "https://files.pythonhosted.org/packages/55/2f/e0e5b237bffdb5d6c530ce87cc1d413a5b7d7dfd60fb067ad6d254c35c76/coverage-7.13.5-cp311-cp311-manylinux1_i686.manylinux_2_28_i686.manylinux_2_5_i686.whl", hash = "sha256:0672854dc733c342fa3e957e0605256d2bf5934feeac328da9e0b5449634a642", size = 250303 },
-    { url = "https://files.pythonhosted.org/packages/92/be/b1afb692be85b947f3401375851484496134c5554e67e822c35f28bf2fbc/coverage-7.13.5-cp311-cp311-manylinux1_x86_64.manylinux_2_28_x86_64.manylinux_2_5_x86_64.whl", hash = "sha256:ec10e2a42b41c923c2209b846126c6582db5e43a33157e9870ba9fb70dc7854b", size = 252218 },
-    { url = "https://files.pythonhosted.org/packages/da/69/2f47bb6fa1b8d1e3e5d0c4be8ccb4313c63d742476a619418f85740d597b/coverage-7.13.5-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:be3d4bbad9d4b037791794ddeedd7d64a56f5933a2c1373e18e9e568b9141686", size = 254326 },
-    { url = "https://files.pythonhosted.org/packages/d5/d0/79db81da58965bd29dabc8f4ad2a2af70611a57cba9d1ec006f072f30a54/coverage-7.13.5-cp311-cp311-manylinux2014_ppc64le.manylinux_2_17_ppc64le.manylinux_2_28_ppc64le.whl", hash = "sha256:4d2afbc5cc54d286bfb54541aa50b64cdb07a718227168c87b9e2fb8f25e1743", size = 256267 },
-    { url = "https://files.pythonhosted.org/packages/e5/32/d0d7cc8168f91ddab44c0ce4806b969df5f5fdfdbb568eaca2dbc2a04936/coverage-7.13.5-cp311-cp311-manylinux_2_31_riscv64.manylinux_2_39_riscv64.whl", hash = "sha256:3ad050321264c49c2fa67bb599100456fc51d004b82534f379d16445da40fb75", size = 250430 },
-    { url = "https://files.pythonhosted.org/packages/4d/06/a055311d891ddbe231cd69fdd20ea4be6e3603ffebddf8704b8ca8e10a3c/coverage-7.13.5-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:7300c8a6d13335b29bb76d7651c66af6bd8658517c43499f110ddc6717bfc209", size = 252017 },
-    { url = "https://files.pythonhosted.org/packages/d6/f6/d0fd2d21e29a657b5f77a2fe7082e1568158340dceb941954f776dce1b7b/coverage-7.13.5-cp311-cp311-musllinux_1_2_i686.whl", hash = "sha256:eb07647a5738b89baab047f14edd18ded523de60f3b30e75c2acc826f79c839a", size = 250080 },
-    { url = "https://files.pythonhosted.org/packages/4e/ab/0d7fb2efc2e9a5eb7ddcc6e722f834a69b454b7e6e5888c3a8567ecffb31/coverage-7.13.5-cp311-cp311-musllinux_1_2_ppc64le.whl", hash = "sha256:9adb6688e3b53adffefd4a52d72cbd8b02602bfb8f74dcd862337182fd4d1a4e", size = 253843 },
-    { url = "https://files.pythonhosted.org/packages/ba/6f/7467b917bbf5408610178f62a49c0ed4377bb16c1657f689cc61470da8ce/coverage-7.13.5-cp311-cp311-musllinux_1_2_riscv64.whl", hash = "sha256:7c8d4bc913dd70b93488d6c496c77f3aff5ea99a07e36a18f865bca55adef8bd", size = 249802 },
-    { url = "https://files.pythonhosted.org/packages/75/2c/1172fb689df92135f5bfbbd69fc83017a76d24ea2e2f3a1154007e2fb9f8/coverage-7.13.5-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:0e3c426ffc4cd952f54ee9ffbdd10345709ecc78a3ecfd796a57236bfad0b9b8", size = 250707 },
-    { url = "https://files.pythonhosted.org/packages/67/21/9ac389377380a07884e3b48ba7a620fcd9dbfaf1d40565facdc6b36ec9ef/coverage-7.13.5-cp311-cp311-win32.whl", hash = "sha256:259b69bb83ad9894c4b25be2528139eecba9a82646ebdda2d9db1ba28424a6bf", size = 221880 },
-    { url = "https://files.pythonhosted.org/packages/af/7f/4cd8a92531253f9d7c1bbecd9fa1b472907fb54446ca768c59b531248dc5/coverage-7.13.5-cp311-cp311-win_amd64.whl", hash = "sha256:258354455f4e86e3e9d0d17571d522e13b4e1e19bf0f8596bcf9476d61e7d8a9", size = 222816 },
-    { url = "https://files.pythonhosted.org/packages/12/a6/1d3f6155fb0010ca68eba7fe48ca6c9da7385058b77a95848710ecf189b1/coverage-7.13.5-cp311-cp311-win_arm64.whl", hash = "sha256:bff95879c33ec8da99fc9b6fe345ddb5be6414b41d6d1ad1c8f188d26f36e028", size = 221483 },
     { url = "https://files.pythonhosted.org/packages/a0/c3/a396306ba7db865bf96fc1fb3b7fd29bcbf3d829df642e77b13555163cd6/coverage-7.13.5-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:460cf0114c5016fa841214ff5564aa4864f11948da9440bc97e21ad1f4ba1e01", size = 219554 },
     { url = "https://files.pythonhosted.org/packages/a6/16/a68a19e5384e93f811dccc51034b1fd0b865841c390e3c931dcc4699e035/coverage-7.13.5-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:0e223ce4b4ed47f065bfb123687686512e37629be25cc63728557ae7db261422", size = 219908 },
     { url = "https://files.pythonhosted.org/packages/29/72/20b917c6793af3a5ceb7fb9c50033f3ec7865f2911a1416b34a7cfa0813b/coverage-7.13.5-cp312-cp312-manylinux1_i686.manylinux_2_28_i686.manylinux_2_5_i686.whl", hash = "sha256:6e3370441f4513c6252bf042b9c36d22491142385049243253c7e48398a15a9f", size = 251419 },
@@ -373,14 +326,6 @@ version = "4.61.0"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/33/f9/0e84d593c0e12244150280a630999835a64f2852276161b62a0f98318de0/fonttools-4.61.0.tar.gz", hash = "sha256:ec520a1f0c7758d7a858a00f090c1745f6cde6a7c5e76fb70ea4044a15f712e7", size = 3561884 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/fd/be/5aa89cdddf2863d8afbdc19eb8ec5d8d35d40eeeb8e6cf52c5ff1c2dbd33/fonttools-4.61.0-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:a32a16951cbf113d38f1dd8551b277b6e06e0f6f776fece0f99f746d739e1be3", size = 2847553 },
-    { url = "https://files.pythonhosted.org/packages/0d/3e/6ff643b07cead1236a534f51291ae2981721cf419135af5b740c002a66dd/fonttools-4.61.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:328a9c227984bebaf69f3ac9062265f8f6acc7ddf2e4e344c63358579af0aa3d", size = 2388298 },
-    { url = "https://files.pythonhosted.org/packages/c3/15/fca8dfbe7b482e6f240b1aad0ed7c6e2e75e7a28efa3d3a03b570617b5e5/fonttools-4.61.0-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:2f0bafc8a3b3749c69cc610e5aa3da832d39c2a37a68f03d18ec9a02ecaac04a", size = 5054133 },
-    { url = "https://files.pythonhosted.org/packages/6a/a2/821c61c691b21fd09e07528a9a499cc2b075ac83ddb644aa16c9875a64bc/fonttools-4.61.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:b5ca59b7417d149cf24e4c1933c9f44b2957424fc03536f132346d5242e0ebe5", size = 5031410 },
-    { url = "https://files.pythonhosted.org/packages/e8/f6/8b16339e93d03c732c8a23edefe3061b17a5f9107ddc47a3215ecd054cac/fonttools-4.61.0-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:df8cbce85cf482eb01f4551edca978c719f099c623277bda8332e5dbe7dba09d", size = 5030005 },
-    { url = "https://files.pythonhosted.org/packages/ac/eb/d4e150427bdaa147755239c931bbce829a88149ade5bfd8a327afe565567/fonttools-4.61.0-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:7fb5b84f48a6a733ca3d7f41aa9551908ccabe8669ffe79586560abcc00a9cfd", size = 5154026 },
-    { url = "https://files.pythonhosted.org/packages/7f/5f/3dd00ce0dba6759943c707b1830af8c0bcf6f8f1a9fe46cb82e7ac2aaa74/fonttools-4.61.0-cp311-cp311-win32.whl", hash = "sha256:787ef9dfd1ea9fe49573c272412ae5f479d78e671981819538143bec65863865", size = 2276035 },
-    { url = "https://files.pythonhosted.org/packages/4e/44/798c472f096ddf12955eddb98f4f7c906e7497695d04ce073ddf7161d134/fonttools-4.61.0-cp311-cp311-win_amd64.whl", hash = "sha256:14fafda386377b6131d9e448af42d0926bad47e038de0e5ba1d58c25d621f028", size = 2327290 },
     { url = "https://files.pythonhosted.org/packages/00/5d/19e5939f773c7cb05480fe2e881d63870b63ee2b4bdb9a77d55b1d36c7b9/fonttools-4.61.0-cp312-cp312-macosx_10_13_universal2.whl", hash = "sha256:e24a1565c4e57111ec7f4915f8981ecbb61adf66a55f378fdc00e206059fcfef", size = 2846930 },
     { url = "https://files.pythonhosted.org/packages/25/b2/0658faf66f705293bd7e739a4f038302d188d424926be9c59bdad945664b/fonttools-4.61.0-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:e2bfacb5351303cae9f072ccf3fc6ecb437a6f359c0606bae4b1ab6715201d87", size = 2383016 },
     { url = "https://files.pythonhosted.org/packages/29/a3/1fa90b95b690f0d7541f48850adc40e9019374d896c1b8148d15012b2458/fonttools-4.61.0-cp312-cp312-manylinux1_x86_64.manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_5_x86_64.whl", hash = "sha256:0bdcf2e29d65c26299cc3d502f4612365e8b90a939f46cd92d037b6cb7bb544a", size = 4949425 },
@@ -495,19 +440,6 @@ version = "1.4.9"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/5c/3c/85844f1b0feb11ee581ac23fe5fce65cd049a200c1446708cc1b7f922875/kiwisolver-1.4.9.tar.gz", hash = "sha256:c3b22c26c6fd6811b0ae8363b95ca8ce4ea3c202d3d0975b2914310ceb1bcc4d", size = 97564 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/6f/ab/c80b0d5a9d8a1a65f4f815f2afff9798b12c3b9f31f1d304dd233dd920e2/kiwisolver-1.4.9-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:eb14a5da6dc7642b0f3a18f13654847cd8b7a2550e2645a5bda677862b03ba16", size = 124167 },
-    { url = "https://files.pythonhosted.org/packages/a0/c0/27fe1a68a39cf62472a300e2879ffc13c0538546c359b86f149cc19f6ac3/kiwisolver-1.4.9-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:39a219e1c81ae3b103643d2aedb90f1ef22650deb266ff12a19e7773f3e5f089", size = 66579 },
-    { url = "https://files.pythonhosted.org/packages/31/a2/a12a503ac1fd4943c50f9822678e8015a790a13b5490354c68afb8489814/kiwisolver-1.4.9-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:2405a7d98604b87f3fc28b1716783534b1b4b8510d8142adca34ee0bc3c87543", size = 65309 },
-    { url = "https://files.pythonhosted.org/packages/66/e1/e533435c0be77c3f64040d68d7a657771194a63c279f55573188161e81ca/kiwisolver-1.4.9-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:dc1ae486f9abcef254b5618dfb4113dd49f94c68e3e027d03cf0143f3f772b61", size = 1435596 },
-    { url = "https://files.pythonhosted.org/packages/67/1e/51b73c7347f9aabdc7215aa79e8b15299097dc2f8e67dee2b095faca9cb0/kiwisolver-1.4.9-cp311-cp311-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:8a1f570ce4d62d718dce3f179ee78dac3b545ac16c0c04bb363b7607a949c0d1", size = 1246548 },
-    { url = "https://files.pythonhosted.org/packages/21/aa/72a1c5d1e430294f2d32adb9542719cfb441b5da368d09d268c7757af46c/kiwisolver-1.4.9-cp311-cp311-manylinux_2_24_ppc64le.manylinux_2_28_ppc64le.whl", hash = "sha256:cb27e7b78d716c591e88e0a09a2139c6577865d7f2e152488c2cc6257f460872", size = 1263618 },
-    { url = "https://files.pythonhosted.org/packages/a3/af/db1509a9e79dbf4c260ce0cfa3903ea8945f6240e9e59d1e4deb731b1a40/kiwisolver-1.4.9-cp311-cp311-manylinux_2_24_s390x.manylinux_2_28_s390x.whl", hash = "sha256:15163165efc2f627eb9687ea5f3a28137217d217ac4024893d753f46bce9de26", size = 1317437 },
-    { url = "https://files.pythonhosted.org/packages/e0/f2/3ea5ee5d52abacdd12013a94130436e19969fa183faa1e7c7fbc89e9a42f/kiwisolver-1.4.9-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:bdee92c56a71d2b24c33a7d4c2856bd6419d017e08caa7802d2963870e315028", size = 2195742 },
-    { url = "https://files.pythonhosted.org/packages/6f/9b/1efdd3013c2d9a2566aa6a337e9923a00590c516add9a1e89a768a3eb2fc/kiwisolver-1.4.9-cp311-cp311-musllinux_1_2_ppc64le.whl", hash = "sha256:412f287c55a6f54b0650bd9b6dce5aceddb95864a1a90c87af16979d37c89771", size = 2290810 },
-    { url = "https://files.pythonhosted.org/packages/fb/e5/cfdc36109ae4e67361f9bc5b41323648cb24a01b9ade18784657e022e65f/kiwisolver-1.4.9-cp311-cp311-musllinux_1_2_s390x.whl", hash = "sha256:2c93f00dcba2eea70af2be5f11a830a742fe6b579a1d4e00f47760ef13be247a", size = 2461579 },
-    { url = "https://files.pythonhosted.org/packages/62/86/b589e5e86c7610842213994cdea5add00960076bef4ae290c5fa68589cac/kiwisolver-1.4.9-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:f117e1a089d9411663a3207ba874f31be9ac8eaa5b533787024dc07aeb74f464", size = 2268071 },
-    { url = "https://files.pythonhosted.org/packages/3b/c6/f8df8509fd1eee6c622febe54384a96cfaf4d43bf2ccec7a0cc17e4715c9/kiwisolver-1.4.9-cp311-cp311-win_amd64.whl", hash = "sha256:be6a04e6c79819c9a8c2373317d19a96048e5a3f90bec587787e86a1153883c2", size = 73840 },
-    { url = "https://files.pythonhosted.org/packages/e2/2d/16e0581daafd147bc11ac53f032a2b45eabac897f42a338d0a13c1e5c436/kiwisolver-1.4.9-cp311-cp311-win_arm64.whl", hash = "sha256:0ae37737256ba2de764ddc12aed4956460277f00c4996d51a197e72f62f5eec7", size = 65159 },
     { url = "https://files.pythonhosted.org/packages/86/c9/13573a747838aeb1c76e3267620daa054f4152444d1f3d1a2324b78255b5/kiwisolver-1.4.9-cp312-cp312-macosx_10_13_universal2.whl", hash = "sha256:ac5a486ac389dddcc5bef4f365b6ae3ffff2c433324fb38dd35e3fab7c957999", size = 123686 },
     { url = "https://files.pythonhosted.org/packages/51/ea/2ecf727927f103ffd1739271ca19c424d0e65ea473fbaeea1c014aea93f6/kiwisolver-1.4.9-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:f2ba92255faa7309d06fe44c3a4a97efe1c8d640c2a79a5ef728b685762a6fd2", size = 66460 },
     { url = "https://files.pythonhosted.org/packages/5b/5a/51f5464373ce2aeb5194508298a508b6f21d3867f499556263c64c621914/kiwisolver-1.4.9-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:4a2899935e724dd1074cb568ce7ac0dce28b2cd6ab539c8e001a8578eb106d14", size = 64952 },
@@ -572,11 +504,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/99/dd/841e9a66c4715477ea0abc78da039832fbb09dac5c35c58dc4c41a407b8a/kiwisolver-1.4.9-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:aedff62918805fb62d43a4aa2ecd4482c380dc76cd31bd7c8878588a61bd0369", size = 2391835 },
     { url = "https://files.pythonhosted.org/packages/0c/28/4b2e5c47a0da96896fdfdb006340ade064afa1e63675d01ea5ac222b6d52/kiwisolver-1.4.9-cp314-cp314t-win_amd64.whl", hash = "sha256:1fa333e8b2ce4d9660f2cda9c0e1b6bafcfb2457a9d259faa82289e73ec24891", size = 79988 },
     { url = "https://files.pythonhosted.org/packages/80/be/3578e8afd18c88cdf9cb4cffde75a96d2be38c5a903f1ed0ceec061bd09e/kiwisolver-1.4.9-cp314-cp314t-win_arm64.whl", hash = "sha256:4a48a2ce79d65d363597ef7b567ce3d14d68783d2b2263d98db3d9477805ba32", size = 70260 },
-    { url = "https://files.pythonhosted.org/packages/a3/0f/36d89194b5a32c054ce93e586d4049b6c2c22887b0eb229c61c68afd3078/kiwisolver-1.4.9-pp311-pypy311_pp73-macosx_10_15_x86_64.whl", hash = "sha256:720e05574713db64c356e86732c0f3c5252818d05f9df320f0ad8380641acea5", size = 60104 },
-    { url = "https://files.pythonhosted.org/packages/52/ba/4ed75f59e4658fd21fe7dde1fee0ac397c678ec3befba3fe6482d987af87/kiwisolver-1.4.9-pp311-pypy311_pp73-macosx_11_0_arm64.whl", hash = "sha256:17680d737d5335b552994a2008fab4c851bcd7de33094a82067ef3a576ff02fa", size = 58592 },
-    { url = "https://files.pythonhosted.org/packages/33/01/a8ea7c5ea32a9b45ceeaee051a04c8ed4320f5add3c51bfa20879b765b70/kiwisolver-1.4.9-pp311-pypy311_pp73-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:85b5352f94e490c028926ea567fc569c52ec79ce131dadb968d3853e809518c2", size = 80281 },
-    { url = "https://files.pythonhosted.org/packages/da/e3/dbd2ecdce306f1d07a1aaf324817ee993aab7aee9db47ceac757deabafbe/kiwisolver-1.4.9-pp311-pypy311_pp73-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:464415881e4801295659462c49461a24fb107c140de781d55518c4b80cb6790f", size = 78009 },
-    { url = "https://files.pythonhosted.org/packages/da/e9/0d4add7873a73e462aeb45c036a2dead2562b825aa46ba326727b3f31016/kiwisolver-1.4.9-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:fb940820c63a9590d31d88b815e7a3aa5915cad3ce735ab45f0c730b39547de1", size = 73929 },
 ]
 
 [[package]]
@@ -594,17 +521,6 @@ version = "3.0.3"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/7e/99/7690b6d4034fffd95959cbe0c02de8deb3098cc577c67bb6a24fe5d7caa7/markupsafe-3.0.3.tar.gz", hash = "sha256:722695808f4b6457b320fdc131280796bdceb04ab50fe1795cd540799ebe1698", size = 80313 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/08/db/fefacb2136439fc8dd20e797950e749aa1f4997ed584c62cfb8ef7c2be0e/markupsafe-3.0.3-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:1cc7ea17a6824959616c525620e387f6dd30fec8cb44f649e31712db02123dad", size = 11631 },
-    { url = "https://files.pythonhosted.org/packages/e1/2e/5898933336b61975ce9dc04decbc0a7f2fee78c30353c5efba7f2d6ff27a/markupsafe-3.0.3-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:4bd4cd07944443f5a265608cc6aab442e4f74dff8088b0dfc8238647b8f6ae9a", size = 12058 },
-    { url = "https://files.pythonhosted.org/packages/1d/09/adf2df3699d87d1d8184038df46a9c80d78c0148492323f4693df54e17bb/markupsafe-3.0.3-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:6b5420a1d9450023228968e7e6a9ce57f65d148ab56d2313fcd589eee96a7a50", size = 24287 },
-    { url = "https://files.pythonhosted.org/packages/30/ac/0273f6fcb5f42e314c6d8cd99effae6a5354604d461b8d392b5ec9530a54/markupsafe-3.0.3-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0bf2a864d67e76e5c9a34dc26ec616a66b9888e25e7b9460e1c76d3293bd9dbf", size = 22940 },
-    { url = "https://files.pythonhosted.org/packages/19/ae/31c1be199ef767124c042c6c3e904da327a2f7f0cd63a0337e1eca2967a8/markupsafe-3.0.3-cp311-cp311-manylinux_2_31_riscv64.manylinux_2_39_riscv64.whl", hash = "sha256:bc51efed119bc9cfdf792cdeaa4d67e8f6fcccab66ed4bfdd6bde3e59bfcbb2f", size = 21887 },
-    { url = "https://files.pythonhosted.org/packages/b2/76/7edcab99d5349a4532a459e1fe64f0b0467a3365056ae550d3bcf3f79e1e/markupsafe-3.0.3-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:068f375c472b3e7acbe2d5318dea141359e6900156b5b2ba06a30b169086b91a", size = 23692 },
-    { url = "https://files.pythonhosted.org/packages/a4/28/6e74cdd26d7514849143d69f0bf2399f929c37dc2b31e6829fd2045b2765/markupsafe-3.0.3-cp311-cp311-musllinux_1_2_riscv64.whl", hash = "sha256:7be7b61bb172e1ed687f1754f8e7484f1c8019780f6f6b0786e76bb01c2ae115", size = 21471 },
-    { url = "https://files.pythonhosted.org/packages/62/7e/a145f36a5c2945673e590850a6f8014318d5577ed7e5920a4b3448e0865d/markupsafe-3.0.3-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:f9e130248f4462aaa8e2552d547f36ddadbeaa573879158d721bbd33dfe4743a", size = 22923 },
-    { url = "https://files.pythonhosted.org/packages/0f/62/d9c46a7f5c9adbeeeda52f5b8d802e1094e9717705a645efc71b0913a0a8/markupsafe-3.0.3-cp311-cp311-win32.whl", hash = "sha256:0db14f5dafddbb6d9208827849fad01f1a2609380add406671a26386cdf15a19", size = 14572 },
-    { url = "https://files.pythonhosted.org/packages/83/8a/4414c03d3f891739326e1783338e48fb49781cc915b2e0ee052aa490d586/markupsafe-3.0.3-cp311-cp311-win_amd64.whl", hash = "sha256:de8a88e63464af587c950061a5e6a67d3632e36df62b986892331d4620a35c01", size = 15077 },
-    { url = "https://files.pythonhosted.org/packages/35/73/893072b42e6862f319b5207adc9ae06070f095b358655f077f69a35601f0/markupsafe-3.0.3-cp311-cp311-win_arm64.whl", hash = "sha256:3b562dd9e9ea93f13d53989d23a7e775fdfd1066c33494ff43f5418bc8c58a5c", size = 13876 },
     { url = "https://files.pythonhosted.org/packages/5a/72/147da192e38635ada20e0a2e1a51cf8823d2119ce8883f7053879c2199b5/markupsafe-3.0.3-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:d53197da72cc091b024dd97249dfc7794d6a56530370992a5e1a08983ad9230e", size = 11615 },
     { url = "https://files.pythonhosted.org/packages/9a/81/7e4e08678a1f98521201c3079f77db69fb552acd56067661f8c2f534a718/markupsafe-3.0.3-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:1872df69a4de6aead3491198eaf13810b565bdbeec3ae2dc8780f14458ec73ce", size = 12020 },
     { url = "https://files.pythonhosted.org/packages/1e/2c/799f4742efc39633a1b54a92eec4082e4f815314869865d876824c257c1e/markupsafe-3.0.3-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:3a7e8ae81ae39e62a41ec302f972ba6ae23a5c5396c8e60113e9066ef893da0d", size = 24332 },
@@ -679,13 +595,6 @@ dependencies = [
 ]
 sdist = { url = "https://files.pythonhosted.org/packages/8a/76/d3c6e3a13fe484ebe7718d14e269c9569c4eb0020a968a327acb3b9a8fe6/matplotlib-3.10.8.tar.gz", hash = "sha256:2299372c19d56bcd35cf05a2738308758d32b9eaed2371898d8f5bd33f084aa3", size = 34806269 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/f8/86/de7e3a1cdcfc941483af70609edc06b83e7c8a0e0dc9ac325200a3f4d220/matplotlib-3.10.8-cp311-cp311-macosx_10_12_x86_64.whl", hash = "sha256:6be43b667360fef5c754dda5d25a32e6307a03c204f3c0fc5468b78fa87b4160", size = 8251215 },
-    { url = "https://files.pythonhosted.org/packages/fd/14/baad3222f424b19ce6ad243c71de1ad9ec6b2e4eb1e458a48fdc6d120401/matplotlib-3.10.8-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:a2b336e2d91a3d7006864e0990c83b216fcdca64b5a6484912902cef87313d78", size = 8139625 },
-    { url = "https://files.pythonhosted.org/packages/8f/a0/7024215e95d456de5883e6732e708d8187d9753a21d32f8ddb3befc0c445/matplotlib-3.10.8-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:efb30e3baaea72ce5928e32bab719ab4770099079d66726a62b11b1ef7273be4", size = 8712614 },
-    { url = "https://files.pythonhosted.org/packages/5a/f4/b8347351da9a5b3f41e26cf547252d861f685c6867d179a7c9d60ad50189/matplotlib-3.10.8-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:d56a1efd5bfd61486c8bc968fa18734464556f0fb8e51690f4ac25d85cbbbbc2", size = 9540997 },
-    { url = "https://files.pythonhosted.org/packages/9e/c0/c7b914e297efe0bc36917bf216b2acb91044b91e930e878ae12981e461e5/matplotlib-3.10.8-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:238b7ce5717600615c895050239ec955d91f321c209dd110db988500558e70d6", size = 9596825 },
-    { url = "https://files.pythonhosted.org/packages/6f/d3/a4bbc01c237ab710a1f22b4da72f4ff6d77eb4c7735ea9811a94ae239067/matplotlib-3.10.8-cp311-cp311-win_amd64.whl", hash = "sha256:18821ace09c763ec93aef5eeff087ee493a24051936d7b9ebcad9662f66501f9", size = 8135090 },
-    { url = "https://files.pythonhosted.org/packages/89/dd/a0b6588f102beab33ca6f5218b31725216577b2a24172f327eaf6417d5c9/matplotlib-3.10.8-cp311-cp311-win_arm64.whl", hash = "sha256:bab485bcf8b1c7d2060b4fcb6fc368a9e6f4cd754c9c2fea281f4be21df394a2", size = 8012377 },
     { url = "https://files.pythonhosted.org/packages/9e/67/f997cdcbb514012eb0d10cd2b4b332667997fb5ebe26b8d41d04962fa0e6/matplotlib-3.10.8-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:64fcc24778ca0404ce0cb7b6b77ae1f4c7231cdd60e6778f999ee05cbd581b9a", size = 8260453 },
     { url = "https://files.pythonhosted.org/packages/7e/65/07d5f5c7f7c994f12c768708bd2e17a4f01a2b0f44a1c9eccad872433e2e/matplotlib-3.10.8-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:b9a5ca4ac220a0cdd1ba6bcba3608547117d30468fefce49bb26f55c1a3d5c58", size = 8148321 },
     { url = "https://files.pythonhosted.org/packages/3e/f3/c5195b1ae57ef85339fd7285dfb603b22c8b4e79114bae5f4f0fcf688677/matplotlib-3.10.8-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:3ab4aabc72de4ff77b3ec33a6d78a68227bf1123465887f9905ba79184a1cc04", size = 8716944 },
@@ -721,9 +630,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/4d/4b/e7beb6bbd49f6bae727a12b270a2654d13c397576d25bd6786e47033300f/matplotlib-3.10.8-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:595ba4d8fe983b88f0eec8c26a241e16d6376fe1979086232f481f8f3f67494c", size = 9614011 },
     { url = "https://files.pythonhosted.org/packages/7c/e6/76f2813d31f032e65f6f797e3f2f6e4aab95b65015924b1c51370395c28a/matplotlib-3.10.8-cp314-cp314t-win_amd64.whl", hash = "sha256:25d380fe8b1dc32cf8f0b1b448470a77afb195438bafdf1d858bfb876f3edf7b", size = 8362801 },
     { url = "https://files.pythonhosted.org/packages/5d/49/d651878698a0b67f23aa28e17f45a6d6dd3d3f933fa29087fa4ce5947b5a/matplotlib-3.10.8-cp314-cp314t-win_arm64.whl", hash = "sha256:113bb52413ea508ce954a02c10ffd0d565f9c3bc7f2eddc27dfe1731e71c7b5f", size = 8192560 },
-    { url = "https://files.pythonhosted.org/packages/04/30/3afaa31c757f34b7725ab9d2ba8b48b5e89c2019c003e7d0ead143aabc5a/matplotlib-3.10.8-pp311-pypy311_pp73-macosx_10_15_x86_64.whl", hash = "sha256:6da7c2ce169267d0d066adcf63758f0604aa6c3eebf67458930f9d9b79ad1db1", size = 8249198 },
-    { url = "https://files.pythonhosted.org/packages/48/2f/6334aec331f57485a642a7c8be03cb286f29111ae71c46c38b363230063c/matplotlib-3.10.8-pp311-pypy311_pp73-macosx_11_0_arm64.whl", hash = "sha256:9153c3292705be9f9c64498a8872118540c3f4123d1a1c840172edf262c8be4a", size = 8136817 },
-    { url = "https://files.pythonhosted.org/packages/73/e4/6d6f14b2a759c622f191b2d67e9075a3f56aaccb3be4bb9bb6890030d0a0/matplotlib-3.10.8-pp311-pypy311_pp73-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:1ae029229a57cd1e8fe542485f27e7ca7b23aa9e8944ddb4985d0bc444f1eca2", size = 8713867 },
 ]
 
 [[package]]
@@ -942,17 +848,6 @@ version = "2.3.5"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/76/65/21b3bc86aac7b8f2862db1e808f1ea22b028e30a225a34a5ede9bf8678f2/numpy-2.3.5.tar.gz", hash = "sha256:784db1dcdab56bf0517743e746dfb0f885fc68d948aba86eeec2cba234bdf1c0", size = 20584950 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/43/77/84dd1d2e34d7e2792a236ba180b5e8fcc1e3e414e761ce0253f63d7f572e/numpy-2.3.5-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:de5672f4a7b200c15a4127042170a694d4df43c992948f5e1af57f0174beed10", size = 17034641 },
-    { url = "https://files.pythonhosted.org/packages/2a/ea/25e26fa5837106cde46ae7d0b667e20f69cbbc0efd64cba8221411ab26ae/numpy-2.3.5-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:acfd89508504a19ed06ef963ad544ec6664518c863436306153e13e94605c218", size = 12528324 },
-    { url = "https://files.pythonhosted.org/packages/4d/1a/e85f0eea4cf03d6a0228f5c0256b53f2df4bc794706e7df019fc622e47f1/numpy-2.3.5-cp311-cp311-macosx_14_0_arm64.whl", hash = "sha256:ffe22d2b05504f786c867c8395de703937f934272eb67586817b46188b4ded6d", size = 5356872 },
-    { url = "https://files.pythonhosted.org/packages/5c/bb/35ef04afd567f4c989c2060cde39211e4ac5357155c1833bcd1166055c61/numpy-2.3.5-cp311-cp311-macosx_14_0_x86_64.whl", hash = "sha256:872a5cf366aec6bb1147336480fef14c9164b154aeb6542327de4970282cd2f5", size = 6893148 },
-    { url = "https://files.pythonhosted.org/packages/f2/2b/05bbeb06e2dff5eab512dfc678b1cc5ee94d8ac5956a0885c64b6b26252b/numpy-2.3.5-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:3095bdb8dd297e5920b010e96134ed91d852d81d490e787beca7e35ae1d89cf7", size = 14557282 },
-    { url = "https://files.pythonhosted.org/packages/65/fb/2b23769462b34398d9326081fad5655198fcf18966fcb1f1e49db44fbf31/numpy-2.3.5-cp311-cp311-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:8cba086a43d54ca804ce711b2a940b16e452807acebe7852ff327f1ecd49b0d4", size = 16897903 },
-    { url = "https://files.pythonhosted.org/packages/ac/14/085f4cf05fc3f1e8aa95e85404e984ffca9b2275a5dc2b1aae18a67538b8/numpy-2.3.5-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:6cf9b429b21df6b99f4dee7a1218b8b7ffbbe7df8764dc0bd60ce8a0708fed1e", size = 16341672 },
-    { url = "https://files.pythonhosted.org/packages/6f/3b/1f73994904142b2aa290449b3bb99772477b5fd94d787093e4f24f5af763/numpy-2.3.5-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:396084a36abdb603546b119d96528c2f6263921c50df3c8fd7cb28873a237748", size = 18838896 },
-    { url = "https://files.pythonhosted.org/packages/cd/b9/cf6649b2124f288309ffc353070792caf42ad69047dcc60da85ee85fea58/numpy-2.3.5-cp311-cp311-win32.whl", hash = "sha256:b0c7088a73aef3d687c4deef8452a3ac7c1be4e29ed8bf3b366c8111128ac60c", size = 6563608 },
-    { url = "https://files.pythonhosted.org/packages/aa/44/9fe81ae1dcc29c531843852e2874080dc441338574ccc4306b39e2ff6e59/numpy-2.3.5-cp311-cp311-win_amd64.whl", hash = "sha256:a414504bef8945eae5f2d7cb7be2d4af77c5d1cb5e20b296c2c25b61dff2900c", size = 13078442 },
-    { url = "https://files.pythonhosted.org/packages/6d/a7/f99a41553d2da82a20a2f22e93c94f928e4490bb447c9ff3c4ff230581d3/numpy-2.3.5-cp311-cp311-win_arm64.whl", hash = "sha256:0cd00b7b36e35398fa2d16af7b907b65304ef8bb4817a550e06e5012929830fa", size = 10458555 },
     { url = "https://files.pythonhosted.org/packages/44/37/e669fe6cbb2b96c62f6bbedc6a81c0f3b7362f6a59230b23caa673a85721/numpy-2.3.5-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:74ae7b798248fe62021dbf3c914245ad45d1a6b0cb4a29ecb4b31d0bfbc4cc3e", size = 16733873 },
     { url = "https://files.pythonhosted.org/packages/c5/65/df0db6c097892c9380851ab9e44b52d4f7ba576b833996e0080181c0c439/numpy-2.3.5-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:ee3888d9ff7c14604052b2ca5535a30216aa0a58e948cdd3eeb8d3415f638769", size = 12259838 },
     { url = "https://files.pythonhosted.org/packages/5b/e1/1ee06e70eb2136797abe847d386e7c0e830b67ad1d43f364dd04fa50d338/numpy-2.3.5-cp312-cp312-macosx_14_0_arm64.whl", hash = "sha256:612a95a17655e213502f60cfb9bf9408efdc9eb1d5f50535cc6eb365d11b42b5", size = 5088378 },
@@ -1008,13 +903,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/51/41/851c4b4082402d9ea860c3626db5d5df47164a712cb23b54be028b184c1c/numpy-2.3.5-cp314-cp314t-win32.whl", hash = "sha256:93eebbcf1aafdf7e2ddd44c2923e2672e1010bddc014138b229e49725b4d6be5", size = 6479806 },
     { url = "https://files.pythonhosted.org/packages/90/30/d48bde1dfd93332fa557cff1972fbc039e055a52021fbef4c2c4b1eefd17/numpy-2.3.5-cp314-cp314t-win_amd64.whl", hash = "sha256:c8a9958e88b65c3b27e22ca2a076311636850b612d6bbfb76e8d156aacde2aaf", size = 13105760 },
     { url = "https://files.pythonhosted.org/packages/2d/fd/4b5eb0b3e888d86aee4d198c23acec7d214baaf17ea93c1adec94c9518b9/numpy-2.3.5-cp314-cp314t-win_arm64.whl", hash = "sha256:6203fdf9f3dc5bdaed7319ad8698e685c7a3be10819f41d32a0723e611733b42", size = 10545459 },
-    { url = "https://files.pythonhosted.org/packages/c6/65/f9dea8e109371ade9c782b4e4756a82edf9d3366bca495d84d79859a0b79/numpy-2.3.5-pp311-pypy311_pp73-macosx_10_15_x86_64.whl", hash = "sha256:f0963b55cdd70fad460fa4c1341f12f976bb26cb66021a5580329bd498988310", size = 16910689 },
-    { url = "https://files.pythonhosted.org/packages/00/4f/edb00032a8fb92ec0a679d3830368355da91a69cab6f3e9c21b64d0bb986/numpy-2.3.5-pp311-pypy311_pp73-macosx_11_0_arm64.whl", hash = "sha256:f4255143f5160d0de972d28c8f9665d882b5f61309d8362fdd3e103cf7bf010c", size = 12457053 },
-    { url = "https://files.pythonhosted.org/packages/16/a4/e8a53b5abd500a63836a29ebe145fc1ab1f2eefe1cfe59276020373ae0aa/numpy-2.3.5-pp311-pypy311_pp73-macosx_14_0_arm64.whl", hash = "sha256:a4b9159734b326535f4dd01d947f919c6eefd2d9827466a696c44ced82dfbc18", size = 5285635 },
-    { url = "https://files.pythonhosted.org/packages/a3/2f/37eeb9014d9c8b3e9c55bc599c68263ca44fdbc12a93e45a21d1d56df737/numpy-2.3.5-pp311-pypy311_pp73-macosx_14_0_x86_64.whl", hash = "sha256:2feae0d2c91d46e59fcd62784a3a83b3fb677fead592ce51b5a6fbb4f95965ff", size = 6801770 },
-    { url = "https://files.pythonhosted.org/packages/7d/e4/68d2f474df2cb671b2b6c2986a02e520671295647dad82484cde80ca427b/numpy-2.3.5-pp311-pypy311_pp73-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:ffac52f28a7849ad7576293c0cb7b9f08304e8f7d738a8cb8a90ec4c55a998eb", size = 14391768 },
-    { url = "https://files.pythonhosted.org/packages/b8/50/94ccd8a2b141cb50651fddd4f6a48874acb3c91c8f0842b08a6afc4b0b21/numpy-2.3.5-pp311-pypy311_pp73-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:63c0e9e7eea69588479ebf4a8a270d5ac22763cc5854e9a7eae952a3908103f7", size = 16729263 },
-    { url = "https://files.pythonhosted.org/packages/2d/ee/346fa473e666fe14c52fcdd19ec2424157290a032d4c41f98127bfb31ac7/numpy-2.3.5-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:f16417ec91f12f814b10bafe79ef77e70113a2f5f7018640e7425ff979253425", size = 12967213 },
 ]
 
 [[package]]
@@ -1046,14 +934,6 @@ dependencies = [
 ]
 sdist = { url = "https://files.pythonhosted.org/packages/da/99/b342345300f13440fe9fe385c3c481e2d9a595ee3bab4d3219247ac94e9a/pandas-3.0.2.tar.gz", hash = "sha256:f4753e73e34c8d83221ba58f232433fca2748be8b18dbca02d242ed153945043", size = 4645855 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/97/35/6411db530c618e0e0005187e35aa02ce60ae4c4c4d206964a2f978217c27/pandas-3.0.2-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:a727a73cbdba2f7458dc82449e2315899d5140b449015d822f515749a46cbbe0", size = 10326926 },
-    { url = "https://files.pythonhosted.org/packages/c4/d3/b7da1d5d7dbdc5ef52ed7debd2b484313b832982266905315dad5a0bf0b1/pandas-3.0.2-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:dbbd4aa20ca51e63b53bbde6a0fa4254b1aaabb74d2f542df7a7959feb1d760c", size = 9926987 },
-    { url = "https://files.pythonhosted.org/packages/52/77/9b1c2d6070b5dbe239a7bc889e21bfa58720793fb902d1e070695d87c6d0/pandas-3.0.2-cp311-cp311-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:339dda302bd8369dedeae979cb750e484d549b563c3f54f3922cb8ff4978c5eb", size = 10757067 },
-    { url = "https://files.pythonhosted.org/packages/20/17/ec40d981705654853726e7ac9aea9ddbb4a5d9cf54d8472222f4f3de06c2/pandas-3.0.2-cp311-cp311-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:61c2fd96d72b983a9891b2598f286befd4ad262161a609c92dc1652544b46b76", size = 11258787 },
-    { url = "https://files.pythonhosted.org/packages/90/e3/3f1126d43d3702ca8773871a81c9f15122a1f412342cc56284ffda5b1f70/pandas-3.0.2-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:c934008c733b8bbea273ea308b73b3156f0181e5b72960790b09c18a2794fe1e", size = 11771616 },
-    { url = "https://files.pythonhosted.org/packages/2e/cf/0f4e268e1f5062e44a6bda9f925806721cd4c95c2b808a4c82ebe914f96b/pandas-3.0.2-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:60a80bb4feacbef5e1447a3f82c33209c8b7e07f28d805cfd1fb951e5cb443aa", size = 12337623 },
-    { url = "https://files.pythonhosted.org/packages/44/a0/97a6339859d4acb2536efb24feb6708e82f7d33b2ed7e036f2983fcced82/pandas-3.0.2-cp311-cp311-win_amd64.whl", hash = "sha256:ed72cb3f45190874eb579c64fa92d9df74e98fd63e2be7f62bce5ace0ade61df", size = 9897372 },
-    { url = "https://files.pythonhosted.org/packages/8f/eb/781516b808a99ddf288143cec46b342b3016c3414d137da1fdc3290d8860/pandas-3.0.2-cp311-cp311-win_arm64.whl", hash = "sha256:f12b1a9e332c01e09510586f8ca9b108fd631fd656af82e452d7315ef6df5f9f", size = 9154922 },
     { url = "https://files.pythonhosted.org/packages/f3/b0/c20bd4d6d3f736e6bd6b55794e9cd0a617b858eaad27c8f410ea05d953b7/pandas-3.0.2-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:232a70ebb568c0c4d2db4584f338c1577d81e3af63292208d615907b698a0f18", size = 10347921 },
     { url = "https://files.pythonhosted.org/packages/35/d0/4831af68ce30cc2d03c697bea8450e3225a835ef497d0d70f31b8cdde965/pandas-3.0.2-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:970762605cff1ca0d3f71ed4f3a769ea8f85fc8e6348f6e110b8fea7e6eb5a14", size = 9888127 },
     { url = "https://files.pythonhosted.org/packages/61/a9/16ea9346e1fc4a96e2896242d9bc674764fb9049b0044c0132502f7a771e/pandas-3.0.2-cp312-cp312-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:aff4e6f4d722e0652707d7bcb190c445fe58428500c6d16005b02401764b1b3d", size = 10399577 },
@@ -1110,17 +990,6 @@ version = "12.0.0"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/5a/b0/cace85a1b0c9775a9f8f5d5423c8261c858760e2466c79b2dd184638b056/pillow-12.0.0.tar.gz", hash = "sha256:87d4f8125c9988bfbed67af47dd7a953e2fc7b0cc1e7800ec6d2080d490bb353", size = 47008828 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/0e/5a/a2f6773b64edb921a756eb0729068acad9fc5208a53f4a349396e9436721/pillow-12.0.0-cp311-cp311-macosx_10_10_x86_64.whl", hash = "sha256:0fd00cac9c03256c8b2ff58f162ebcd2587ad3e1f2e397eab718c47e24d231cc", size = 5289798 },
-    { url = "https://files.pythonhosted.org/packages/2e/05/069b1f8a2e4b5a37493da6c5868531c3f77b85e716ad7a590ef87d58730d/pillow-12.0.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:a3475b96f5908b3b16c47533daaa87380c491357d197564e0ba34ae75c0f3257", size = 4650589 },
-    { url = "https://files.pythonhosted.org/packages/61/e3/2c820d6e9a36432503ead175ae294f96861b07600a7156154a086ba7111a/pillow-12.0.0-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:110486b79f2d112cf6add83b28b627e369219388f64ef2f960fef9ebaf54c642", size = 6230472 },
-    { url = "https://files.pythonhosted.org/packages/4f/89/63427f51c64209c5e23d4d52071c8d0f21024d3a8a487737caaf614a5795/pillow-12.0.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:5269cc1caeedb67e6f7269a42014f381f45e2e7cd42d834ede3c703a1d915fe3", size = 8033887 },
-    { url = "https://files.pythonhosted.org/packages/f6/1b/c9711318d4901093c15840f268ad649459cd81984c9ec9887756cca049a5/pillow-12.0.0-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:aa5129de4e174daccbc59d0a3b6d20eaf24417d59851c07ebb37aeb02947987c", size = 6343964 },
-    { url = "https://files.pythonhosted.org/packages/41/1e/db9470f2d030b4995083044cd8738cdd1bf773106819f6d8ba12597d5352/pillow-12.0.0-cp311-cp311-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:bee2a6db3a7242ea309aa7ee8e2780726fed67ff4e5b40169f2c940e7eb09227", size = 7034756 },
-    { url = "https://files.pythonhosted.org/packages/cc/b0/6177a8bdd5ee4ed87cba2de5a3cc1db55ffbbec6176784ce5bb75aa96798/pillow-12.0.0-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:90387104ee8400a7b4598253b4c406f8958f59fcf983a6cea2b50d59f7d63d0b", size = 6458075 },
-    { url = "https://files.pythonhosted.org/packages/bc/5e/61537aa6fa977922c6a03253a0e727e6e4a72381a80d63ad8eec350684f2/pillow-12.0.0-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:bc91a56697869546d1b8f0a3ff35224557ae7f881050e99f615e0119bf934b4e", size = 7125955 },
-    { url = "https://files.pythonhosted.org/packages/1f/3d/d5033539344ee3cbd9a4d69e12e63ca3a44a739eb2d4c8da350a3d38edd7/pillow-12.0.0-cp311-cp311-win32.whl", hash = "sha256:27f95b12453d165099c84f8a8bfdfd46b9e4bda9e0e4b65f0635430027f55739", size = 6298440 },
-    { url = "https://files.pythonhosted.org/packages/4d/42/aaca386de5cc8bd8a0254516957c1f265e3521c91515b16e286c662854c4/pillow-12.0.0-cp311-cp311-win_amd64.whl", hash = "sha256:b583dc9070312190192631373c6c8ed277254aa6e6084b74bdd0a6d3b221608e", size = 6999256 },
-    { url = "https://files.pythonhosted.org/packages/ba/f1/9197c9c2d5708b785f631a6dfbfa8eb3fb9672837cb92ae9af812c13b4ed/pillow-12.0.0-cp311-cp311-win_arm64.whl", hash = "sha256:759de84a33be3b178a64c8ba28ad5c135900359e85fb662bc6e403ad4407791d", size = 2436025 },
     { url = "https://files.pythonhosted.org/packages/2c/90/4fcce2c22caf044e660a198d740e7fbc14395619e3cb1abad12192c0826c/pillow-12.0.0-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:53561a4ddc36facb432fae7a9d8afbfaf94795414f5cdc5fc52f28c1dca90371", size = 5249377 },
     { url = "https://files.pythonhosted.org/packages/fd/e0/ed960067543d080691d47d6938ebccbf3976a931c9567ab2fbfab983a5dd/pillow-12.0.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:71db6b4c1653045dacc1585c1b0d184004f0d7e694c7b34ac165ca70c0838082", size = 4650343 },
     { url = "https://files.pythonhosted.org/packages/e7/a1/f81fdeddcb99c044bf7d6faa47e12850f13cee0849537a7d27eeab5534d4/pillow-12.0.0-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:2fa5f0b6716fc88f11380b88b31fe591a06c6315e955c096c35715788b339e3f", size = 6232981 },
@@ -1182,13 +1051,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/a2/2f/16cabcc6426c32218ace36bf0d55955e813f2958afddbf1d391849fee9d1/pillow-12.0.0-cp314-cp314t-win32.whl", hash = "sha256:3830c769decf88f1289680a59d4f4c46c72573446352e2befec9a8512104fa52", size = 6408045 },
     { url = "https://files.pythonhosted.org/packages/35/73/e29aa0c9c666cf787628d3f0dcf379f4791fba79f4936d02f8b37165bdf8/pillow-12.0.0-cp314-cp314t-win_amd64.whl", hash = "sha256:905b0365b210c73afb0ebe9101a32572152dfd1c144c7e28968a331b9217b94a", size = 7148282 },
     { url = "https://files.pythonhosted.org/packages/c1/70/6b41bdcddf541b437bbb9f47f94d2db5d9ddef6c37ccab8c9107743748a4/pillow-12.0.0-cp314-cp314t-win_arm64.whl", hash = "sha256:99353a06902c2e43b43e8ff74ee65a7d90307d82370604746738a1e0661ccca7", size = 2525630 },
-    { url = "https://files.pythonhosted.org/packages/1d/b3/582327e6c9f86d037b63beebe981425d6811104cb443e8193824ef1a2f27/pillow-12.0.0-pp311-pypy311_pp73-macosx_10_15_x86_64.whl", hash = "sha256:b22bd8c974942477156be55a768f7aa37c46904c175be4e158b6a86e3a6b7ca8", size = 5215068 },
-    { url = "https://files.pythonhosted.org/packages/fd/d6/67748211d119f3b6540baf90f92fae73ae51d5217b171b0e8b5f7e5d558f/pillow-12.0.0-pp311-pypy311_pp73-macosx_11_0_arm64.whl", hash = "sha256:805ebf596939e48dbb2e4922a1d3852cfc25c38160751ce02da93058b48d252a", size = 4614994 },
-    { url = "https://files.pythonhosted.org/packages/2d/e1/f8281e5d844c41872b273b9f2c34a4bf64ca08905668c8ae730eedc7c9fa/pillow-12.0.0-pp311-pypy311_pp73-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:cae81479f77420d217def5f54b5b9d279804d17e982e0f2fa19b1d1e14ab5197", size = 5246639 },
-    { url = "https://files.pythonhosted.org/packages/94/5a/0d8ab8ffe8a102ff5df60d0de5af309015163bf710c7bb3e8311dd3b3ad0/pillow-12.0.0-pp311-pypy311_pp73-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:aeaefa96c768fc66818730b952a862235d68825c178f1b3ffd4efd7ad2edcb7c", size = 6986839 },
-    { url = "https://files.pythonhosted.org/packages/20/2e/3434380e8110b76cd9eb00a363c484b050f949b4bbe84ba770bb8508a02c/pillow-12.0.0-pp311-pypy311_pp73-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:09f2d0abef9e4e2f349305a4f8cc784a8a6c2f58a8c4892eea13b10a943bd26e", size = 5313505 },
-    { url = "https://files.pythonhosted.org/packages/57/ca/5a9d38900d9d74785141d6580950fe705de68af735ff6e727cb911b64740/pillow-12.0.0-pp311-pypy311_pp73-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:bdee52571a343d721fb2eb3b090a82d959ff37fc631e3f70422e0c2e029f3e76", size = 5963654 },
-    { url = "https://files.pythonhosted.org/packages/95/7e/f896623c3c635a90537ac093c6a618ebe1a90d87206e42309cb5d98a1b9e/pillow-12.0.0-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:b290fd8aa38422444d4b50d579de197557f182ef1068b75f5aa8558638b8d0a5", size = 6997850 },
 ]
 
 [[package]]
@@ -1274,15 +1136,6 @@ version = "6.0.3"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/05/8e/961c0007c59b8dd7729d542c61a4d537767a59645b82a0b521206e1e25c2/pyyaml-6.0.3.tar.gz", hash = "sha256:d76623373421df22fb4cf8817020cbb7ef15c725b9d5e45f17e189bfc384190f", size = 130960 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/6d/16/a95b6757765b7b031c9374925bb718d55e0a9ba8a1b6a12d25962ea44347/pyyaml-6.0.3-cp311-cp311-macosx_10_13_x86_64.whl", hash = "sha256:44edc647873928551a01e7a563d7452ccdebee747728c1080d881d68af7b997e", size = 185826 },
-    { url = "https://files.pythonhosted.org/packages/16/19/13de8e4377ed53079ee996e1ab0a9c33ec2faf808a4647b7b4c0d46dd239/pyyaml-6.0.3-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:652cb6edd41e718550aad172851962662ff2681490a8a711af6a4d288dd96824", size = 175577 },
-    { url = "https://files.pythonhosted.org/packages/0c/62/d2eb46264d4b157dae1275b573017abec435397aa59cbcdab6fc978a8af4/pyyaml-6.0.3-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:10892704fc220243f5305762e276552a0395f7beb4dbf9b14ec8fd43b57f126c", size = 775556 },
-    { url = "https://files.pythonhosted.org/packages/10/cb/16c3f2cf3266edd25aaa00d6c4350381c8b012ed6f5276675b9eba8d9ff4/pyyaml-6.0.3-cp311-cp311-manylinux2014_s390x.manylinux_2_17_s390x.manylinux_2_28_s390x.whl", hash = "sha256:850774a7879607d3a6f50d36d04f00ee69e7fc816450e5f7e58d7f17f1ae5c00", size = 882114 },
-    { url = "https://files.pythonhosted.org/packages/71/60/917329f640924b18ff085ab889a11c763e0b573da888e8404ff486657602/pyyaml-6.0.3-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:b8bb0864c5a28024fac8a632c443c87c5aa6f215c0b126c449ae1a150412f31d", size = 806638 },
-    { url = "https://files.pythonhosted.org/packages/dd/6f/529b0f316a9fd167281a6c3826b5583e6192dba792dd55e3203d3f8e655a/pyyaml-6.0.3-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:1d37d57ad971609cf3c53ba6a7e365e40660e3be0e5175fa9f2365a379d6095a", size = 767463 },
-    { url = "https://files.pythonhosted.org/packages/f2/6a/b627b4e0c1dd03718543519ffb2f1deea4a1e6d42fbab8021936a4d22589/pyyaml-6.0.3-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:37503bfbfc9d2c40b344d06b2199cf0e96e97957ab1c1b546fd4f87e53e5d3e4", size = 794986 },
-    { url = "https://files.pythonhosted.org/packages/45/91/47a6e1c42d9ee337c4839208f30d9f09caa9f720ec7582917b264defc875/pyyaml-6.0.3-cp311-cp311-win32.whl", hash = "sha256:8098f252adfa6c80ab48096053f512f2321f0b998f98150cea9bd23d83e1467b", size = 142543 },
-    { url = "https://files.pythonhosted.org/packages/da/e3/ea007450a105ae919a72393cb06f122f288ef60bba2dc64b26e2646fa315/pyyaml-6.0.3-cp311-cp311-win_amd64.whl", hash = "sha256:9f3bfb4965eb874431221a3ff3fdcddc7e74e3b07799e0e84ca4a0f867d449bf", size = 158763 },
     { url = "https://files.pythonhosted.org/packages/d1/33/422b98d2195232ca1826284a76852ad5a86fe23e31b009c9886b2d0fb8b2/pyyaml-6.0.3-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:7f047e29dcae44602496db43be01ad42fc6f1cc0d8cd6c83d342306c32270196", size = 182063 },
     { url = "https://files.pythonhosted.org/packages/89/a0/6cf41a19a1f2f3feab0e9c0b74134aa2ce6849093d5517a0c550fe37a648/pyyaml-6.0.3-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:fc09d0aa354569bc501d4e787133afc08552722d3ab34836a80547331bb5d4a0", size = 173973 },
     { url = "https://files.pythonhosted.org/packages/ed/23/7a778b6bd0b9a8039df8b1b1d80e2e2ad78aa04171592c8a5c43a56a6af4/pyyaml-6.0.3-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:9149cad251584d5fb4981be1ecde53a1ca46c891a79788c0df828d2f166bda28", size = 775116 },
@@ -1399,16 +1252,6 @@ dependencies = [
 ]
 sdist = { url = "https://files.pythonhosted.org/packages/0a/ca/d8ace4f98322d01abcd52d381134344bf7b431eba7ed8b42bdea5a3c2ac9/scipy-1.16.3.tar.gz", hash = "sha256:01e87659402762f43bd2fee13370553a17ada367d42e7487800bf2916535aecb", size = 30597883 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/9b/5f/6f37d7439de1455ce9c5a556b8d1db0979f03a796c030bafdf08d35b7bf9/scipy-1.16.3-cp311-cp311-macosx_10_14_x86_64.whl", hash = "sha256:40be6cf99e68b6c4321e9f8782e7d5ff8265af28ef2cd56e9c9b2638fa08ad97", size = 36630881 },
-    { url = "https://files.pythonhosted.org/packages/7c/89/d70e9f628749b7e4db2aa4cd89735502ff3f08f7b9b27d2e799485987cd9/scipy-1.16.3-cp311-cp311-macosx_12_0_arm64.whl", hash = "sha256:8be1ca9170fcb6223cc7c27f4305d680ded114a1567c0bd2bfcbf947d1b17511", size = 28941012 },
-    { url = "https://files.pythonhosted.org/packages/a8/a8/0e7a9a6872a923505dbdf6bb93451edcac120363131c19013044a1e7cb0c/scipy-1.16.3-cp311-cp311-macosx_14_0_arm64.whl", hash = "sha256:bea0a62734d20d67608660f69dcda23e7f90fb4ca20974ab80b6ed40df87a005", size = 20931935 },
-    { url = "https://files.pythonhosted.org/packages/bd/c7/020fb72bd79ad798e4dbe53938543ecb96b3a9ac3fe274b7189e23e27353/scipy-1.16.3-cp311-cp311-macosx_14_0_x86_64.whl", hash = "sha256:2a207a6ce9c24f1951241f4693ede2d393f59c07abc159b2cb2be980820e01fb", size = 23534466 },
-    { url = "https://files.pythonhosted.org/packages/be/a0/668c4609ce6dbf2f948e167836ccaf897f95fb63fa231c87da7558a374cd/scipy-1.16.3-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:532fb5ad6a87e9e9cd9c959b106b73145a03f04c7d57ea3e6f6bb60b86ab0876", size = 33593618 },
-    { url = "https://files.pythonhosted.org/packages/ca/6e/8942461cf2636cdae083e3eb72622a7fbbfa5cf559c7d13ab250a5dbdc01/scipy-1.16.3-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:0151a0749efeaaab78711c78422d413c583b8cdd2011a3c1d6c794938ee9fdb2", size = 35899798 },
-    { url = "https://files.pythonhosted.org/packages/79/e8/d0f33590364cdbd67f28ce79368b373889faa4ee959588beddf6daef9abe/scipy-1.16.3-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:b7180967113560cca57418a7bc719e30366b47959dd845a93206fbed693c867e", size = 36226154 },
-    { url = "https://files.pythonhosted.org/packages/39/c1/1903de608c0c924a1749c590064e65810f8046e437aba6be365abc4f7557/scipy-1.16.3-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:deb3841c925eeddb6afc1e4e4a45e418d19ec7b87c5df177695224078e8ec733", size = 38878540 },
-    { url = "https://files.pythonhosted.org/packages/f1/d0/22ec7036ba0b0a35bccb7f25ab407382ed34af0b111475eb301c16f8a2e5/scipy-1.16.3-cp311-cp311-win_amd64.whl", hash = "sha256:53c3844d527213631e886621df5695d35e4f6a75f620dca412bcd292f6b87d78", size = 38722107 },
-    { url = "https://files.pythonhosted.org/packages/7b/60/8a00e5a524bb3bf8898db1650d350f50e6cffb9d7a491c561dc9826c7515/scipy-1.16.3-cp311-cp311-win_arm64.whl", hash = "sha256:9452781bd879b14b6f055b26643703551320aa8d79ae064a71df55c00286a184", size = 25506272 },
     { url = "https://files.pythonhosted.org/packages/40/41/5bf55c3f386b1643812f3a5674edf74b26184378ef0f3e7c7a09a7e2ca7f/scipy-1.16.3-cp312-cp312-macosx_10_14_x86_64.whl", hash = "sha256:81fc5827606858cf71446a5e98715ba0e11f0dbc83d71c7409d05486592a45d6", size = 36659043 },
     { url = "https://files.pythonhosted.org/packages/1e/0f/65582071948cfc45d43e9870bf7ca5f0e0684e165d7c9ef4e50d783073eb/scipy-1.16.3-cp312-cp312-macosx_12_0_arm64.whl", hash = "sha256:c97176013d404c7346bf57874eaac5187d969293bf40497140b0a2b2b7482e07", size = 28898986 },
     { url = "https://files.pythonhosted.org/packages/96/5e/36bf3f0ac298187d1ceadde9051177d6a4fe4d507e8f59067dc9dd39e650/scipy-1.16.3-cp312-cp312-macosx_14_0_arm64.whl", hash = "sha256:2b71d93c8a9936046866acebc915e2af2e292b883ed6e2cbe5c34beb094b82d9", size = 20889814 },
@@ -1512,9 +1355,6 @@ version = "6.0.0"
 source = { registry = "https://pypi.org/simple" }
 sdist = { url = "https://files.pythonhosted.org/packages/db/7d/7f3d619e951c88ed75c6037b246ddcf2d322812ee8ea189be89511721d54/watchdog-6.0.0.tar.gz", hash = "sha256:9ddf7c82fda3ae8e24decda1338ede66e1c99883db93711d8fb941eaa2d8c282", size = 131220 }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/e0/24/d9be5cd6642a6aa68352ded4b4b10fb0d7889cb7f45814fb92cecd35f101/watchdog-6.0.0-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:6eb11feb5a0d452ee41f824e271ca311a09e250441c262ca2fd7ebcf2461a06c", size = 96393 },
-    { url = "https://files.pythonhosted.org/packages/63/7a/6013b0d8dbc56adca7fdd4f0beed381c59f6752341b12fa0886fa7afc78b/watchdog-6.0.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:ef810fbf7b781a5a593894e4f439773830bdecb885e6880d957d5b9382a960d2", size = 88392 },
-    { url = "https://files.pythonhosted.org/packages/d1/40/b75381494851556de56281e053700e46bff5b37bf4c7267e858640af5a7f/watchdog-6.0.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:afd0fe1b2270917c5e23c2a65ce50c2a4abb63daafb0d419fde368e272a76b7c", size = 89019 },
     { url = "https://files.pythonhosted.org/packages/39/ea/3930d07dafc9e286ed356a679aa02d777c06e9bfd1164fa7c19c288a5483/watchdog-6.0.0-cp312-cp312-macosx_10_13_universal2.whl", hash = "sha256:bdd4e6f14b8b18c334febb9c4425a878a2ac20efd1e0b231978e7b150f92a948", size = 96471 },
     { url = "https://files.pythonhosted.org/packages/12/87/48361531f70b1f87928b045df868a9fd4e253d9ae087fa4cf3f7113be363/watchdog-6.0.0-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:c7c15dda13c4eb00d6fb6fc508b3c0ed88b9d5d374056b239c4ad1611125c860", size = 88449 },
     { url = "https://files.pythonhosted.org/packages/5b/7e/8f322f5e600812e6f9a31b75d242631068ca8f4ef0582dd3ae6e72daecc8/watchdog-6.0.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:6f10cb2d5902447c7d0da897e2c6768bca89174d0c6e1e30abec5421af97a5b0", size = 89054 },