From 267d161769e709a7631b1134efbf08b1ad2908c0 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Fri, 17 Apr 2026 23:25:38 -0700 Subject: [PATCH] [tests] more test coverage --- meanas/fdtd/pml.py | 4 - meanas/test/test_bloch_interactions.py | 189 +++++++++++++++++++++ meanas/test/test_eigensolvers_numerics.py | 36 ++++ meanas/test/test_fdfd_algebra_helpers.py | 9 + meanas/test/test_fdfd_farfield.py | 74 ++++++++ meanas/test/test_fdtd_base.py | 44 +++++ meanas/test/test_import_fallbacks.py | 43 +++++ meanas/test/test_waveguide_mode_helpers.py | 15 ++ 8 files changed, 410 insertions(+), 4 deletions(-) create mode 100644 meanas/test/test_fdfd_farfield.py create mode 100644 meanas/test/test_fdtd_base.py create mode 100644 meanas/test/test_import_fallbacks.py diff --git a/meanas/fdtd/pml.py b/meanas/fdtd/pml.py index fc456eb..dec3d83 100644 --- a/meanas/fdtd/pml.py +++ b/meanas/fdtd/pml.py @@ -57,8 +57,6 @@ def cpml_params( xh -= 0.5 xe = xe[::-1] xh = xh[::-1] - else: - raise Exception('Bad polarity!') expand_slice_l: list[Any] = [None, None, None] expand_slice_l[axis] = slice(None) @@ -82,8 +80,6 @@ def cpml_params( region_list[axis] = slice(None, thickness) elif polarity > 0: region_list[axis] = slice(-thickness, None) - else: - raise Exception('Bad polarity!') region = tuple(region_list) return { diff --git a/meanas/test/test_bloch_interactions.py b/meanas/test/test_bloch_interactions.py index 28edcf8..0a2b70c 100644 --- a/meanas/test/test_bloch_interactions.py +++ b/meanas/test/test_bloch_interactions.py @@ -1,5 +1,7 @@ import numpy +import pytest from numpy.testing import assert_allclose +from types import SimpleNamespace from ..fdfd import bloch @@ -10,6 +12,12 @@ EPSILON = numpy.ones((3, *SHAPE), dtype=float) K0 = numpy.array([0.1, 0.0, 0.0], dtype=float) H_SIZE = 2 * numpy.prod(SHAPE) Y0 = (numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE))[None, :] +Y0_TWO_MODE = numpy.vstack( + [ + numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE), + numpy.linspace(2.0, 3.5, H_SIZE) - 0.5j * numpy.arange(H_SIZE, dtype=float), + ], +) def build_overlap_fixture() -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]: @@ -98,6 +106,187 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None: assert callback_count > 0 +def test_eigsolve_without_initial_guess_returns_finite_modes() -> None: + eigvals, eigvecs = bloch.eigsolve( + 1, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=20, + y0=None, + ) + + operator = bloch.maxwell_operator(K0, G_MATRIX, EPSILON) + eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0]) + residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec) + + assert eigvals.shape == (1,) + assert eigvecs.shape == (1, H_SIZE) + assert numpy.isfinite(eigvals).all() + assert numpy.isfinite(eigvecs).all() + assert residual < 1e-5 + + +def test_eigsolve_recovers_from_singular_initial_subspace(monkeypatch: pytest.MonkeyPatch) -> None: + class FakeRng: + def __init__(self) -> None: + self.calls = 0 + + def random(self, shape: tuple[int, ...]) -> numpy.ndarray: + self.calls += 1 + return numpy.arange(numpy.prod(shape), dtype=float).reshape(shape) + self.calls + + fake_rng = FakeRng() + singular_y0 = numpy.vstack([Y0_TWO_MODE[0], Y0_TWO_MODE[0]]) + monkeypatch.setattr(bloch.numpy.random, 'default_rng', lambda: fake_rng) + + eigvals, eigvecs = bloch.eigsolve( + 2, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=20, + y0=singular_y0, + ) + + assert fake_rng.calls == 2 + assert eigvals.shape == (2,) + assert eigvecs.shape == (2, H_SIZE) + assert numpy.isfinite(eigvals).all() + assert numpy.isfinite(eigvecs).all() + + +def test_eigsolve_reconditions_large_trace_initial_subspace(monkeypatch: pytest.MonkeyPatch) -> None: + original_inv = bloch.numpy.linalg.inv + original_sqrtm = bloch.scipy.linalg.sqrtm + sqrtm_calls = 0 + inv_calls = 0 + + def inv_with_large_first_trace(matrix: numpy.ndarray) -> numpy.ndarray: + nonlocal inv_calls + inv_calls += 1 + if inv_calls == 1: + return numpy.eye(matrix.shape[0], dtype=complex) * 1e9 + return original_inv(matrix) + + def sqrtm_wrapper(matrix: numpy.ndarray) -> numpy.ndarray: + nonlocal sqrtm_calls + sqrtm_calls += 1 + return original_sqrtm(matrix) + + monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_with_large_first_trace) + monkeypatch.setattr(bloch.scipy.linalg, 'sqrtm', sqrtm_wrapper) + + eigvals, eigvecs = bloch.eigsolve( + 2, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=20, + y0=Y0_TWO_MODE, + ) + + assert sqrtm_calls >= 2 + assert eigvals.shape == (2,) + assert eigvecs.shape == (2, H_SIZE) + assert numpy.isfinite(eigvals).all() + assert numpy.isfinite(eigvecs).all() + + +def test_eigsolve_qi_memoization_reuses_cached_theta(monkeypatch: pytest.MonkeyPatch) -> None: + def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace: + theta = 0.3 + first = func(theta) + second = func(theta) + assert_allclose(second, first) + return SimpleNamespace(fun=second, x=theta) + + monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar) + + eigvals, eigvecs = bloch.eigsolve( + 1, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=1, + y0=Y0, + ) + + assert eigvals.shape == (1,) + assert eigvecs.shape == (1, H_SIZE) + assert numpy.isfinite(eigvals).all() + assert numpy.isfinite(eigvecs).all() + + +@pytest.mark.parametrize('theta', [numpy.pi / 2 - 1e-8, 1e-8]) +def test_eigsolve_qi_taylor_expansions_return_finite_modes(monkeypatch: pytest.MonkeyPatch, theta: float) -> None: + original_inv = bloch.numpy.linalg.inv + inv_calls = 0 + + def inv_raise_once_for_q(matrix: numpy.ndarray) -> numpy.ndarray: + nonlocal inv_calls + inv_calls += 1 + if inv_calls == 3: + raise numpy.linalg.LinAlgError('forced singular Q') + return original_inv(matrix) + + def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace: + value = func(theta) + return SimpleNamespace(fun=value, x=theta) + + monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_raise_once_for_q) + monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar) + + eigvals, eigvecs = bloch.eigsolve( + 1, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=1, + y0=Y0, + ) + + assert eigvals.shape == (1,) + assert eigvecs.shape == (1, H_SIZE) + assert numpy.isfinite(eigvals).all() + assert numpy.isfinite(eigvecs).all() + + +def test_eigsolve_qi_inexplicable_singularity_raises(monkeypatch: pytest.MonkeyPatch) -> None: + original_inv = bloch.numpy.linalg.inv + inv_calls = 0 + + def inv_raise_once_for_q(matrix: numpy.ndarray) -> numpy.ndarray: + nonlocal inv_calls + inv_calls += 1 + if inv_calls == 3: + raise numpy.linalg.LinAlgError('forced singular Q') + return original_inv(matrix) + + def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace: + func(numpy.pi / 4) + raise AssertionError('unreachable after trace_func exception') + + monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_raise_once_for_q) + monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar) + + with pytest.raises(Exception, match='Inexplicable singularity in trace_func'): + bloch.eigsolve( + 1, + K0, + G_MATRIX, + EPSILON, + tolerance=1e-6, + max_iters=1, + y0=Y0, + ) + + def test_find_k_returns_vector_frequency_and_callbacks() -> None: target_eigvals, _target_eigvecs = bloch.eigsolve( 1, diff --git a/meanas/test/test_eigensolvers_numerics.py b/meanas/test/test_eigensolvers_numerics.py index 8d0c5c9..5849d2c 100644 --- a/meanas/test/test_eigensolvers_numerics.py +++ b/meanas/test/test_eigensolvers_numerics.py @@ -2,6 +2,7 @@ import numpy from numpy.linalg import norm from scipy import sparse import scipy.sparse.linalg as spalg +from numpy.testing import assert_allclose from ..eigensolvers import rayleigh_quotient_iteration, signed_eigensolve @@ -45,3 +46,38 @@ def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None: assert eigvecs.shape == (4, 1) assert abs(eigvals[0] + 2.0) < 1e-12 assert abs(eigvecs[3, 0]) > 0.99 + + +def test_rayleigh_quotient_iteration_uses_default_linear_operator_solver() -> None: + operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() + linear_operator = spalg.LinearOperator( + shape=operator.shape, + dtype=complex, + matvec=lambda vv: operator @ vv, + ) + + eigval, eigvec = rayleigh_quotient_iteration( + linear_operator, + numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex), + iterations=8, + ) + + residual = norm(operator @ eigvec - eigval * eigvec) + assert abs(eigval - 3.0) < 1e-12 + assert residual < 1e-12 + + +def test_signed_eigensolve_linear_operator_fallback_returns_dominant_positive_mode() -> None: + operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() + linear_operator = spalg.LinearOperator( + shape=operator.shape, + dtype=complex, + matvec=lambda vv: operator @ vv, + ) + + eigvals, eigvecs = signed_eigensolve(linear_operator, how_many=1) + + assert eigvals.shape == (1,) + assert eigvecs.shape == (4, 1) + assert_allclose(eigvals[0], 5.0, atol=1e-12, rtol=1e-12) + assert abs(eigvecs[0, 0]) > 0.99 diff --git a/meanas/test/test_fdfd_algebra_helpers.py b/meanas/test/test_fdfd_algebra_helpers.py index b481023..57bb431 100644 --- a/meanas/test/test_fdfd_algebra_helpers.py +++ b/meanas/test/test_fdfd_algebra_helpers.py @@ -200,6 +200,15 @@ def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None: assert numpy.isfinite(dxes[1][0]).all() +def test_uniform_grid_scpml_default_s_function_matches_explicit_default() -> None: + implicit = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0) + explicit = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0, s_function=scpml.prepare_s_function()) + + for implicit_group, explicit_group in zip(implicit, explicit, strict=True): + for implicit_axis, explicit_axis in zip(implicit_group, explicit_group, strict=True): + assert_allclose(implicit_axis, explicit_axis) + + 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)] diff --git a/meanas/test/test_fdfd_farfield.py b/meanas/test/test_fdfd_farfield.py new file mode 100644 index 0000000..a040b67 --- /dev/null +++ b/meanas/test/test_fdfd_farfield.py @@ -0,0 +1,74 @@ +import numpy +import pytest + +from ..fdfd import farfield + + +NEAR_SHAPE = (2, 3) +E_NEAR = [numpy.zeros(NEAR_SHAPE, dtype=complex), numpy.zeros(NEAR_SHAPE, dtype=complex)] +H_NEAR = [numpy.zeros(NEAR_SHAPE, dtype=complex), numpy.zeros(NEAR_SHAPE, dtype=complex)] + + +def test_near_to_farfield_rejects_wrong_length_inputs() -> None: + with pytest.raises(Exception, match='E_near must be a length-2 list'): + farfield.near_to_farfield(E_NEAR[:1], H_NEAR, dx=0.2, dy=0.3) + + with pytest.raises(Exception, match='H_near must be a length-2 list'): + farfield.near_to_farfield(E_NEAR, H_NEAR[:1], dx=0.2, dy=0.3) + + +def test_near_to_farfield_rejects_mismatched_shapes() -> None: + bad_h_near = [H_NEAR[0], numpy.zeros((2, 4), dtype=complex)] + + with pytest.raises(Exception, match='All fields must be the same shape'): + farfield.near_to_farfield(E_NEAR, bad_h_near, dx=0.2, dy=0.3) + + +def test_near_to_farfield_uses_default_and_scalar_padding_shapes() -> None: + default_result = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3) + scalar_result = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8) + + assert default_result['E'][0].shape == (2, 4) + assert default_result['H'][0].shape == (2, 4) + assert scalar_result['E'][0].shape == (8, 8) + assert scalar_result['H'][0].shape == (8, 8) + + +def test_far_to_nearfield_rejects_wrong_length_inputs() -> None: + ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8) + + with pytest.raises(Exception, match='E_far must be a length-2 list'): + farfield.far_to_nearfield(ff['E'][:1], ff['H'], ff['dkx'], ff['dky']) + + with pytest.raises(Exception, match='H_far must be a length-2 list'): + farfield.far_to_nearfield(ff['E'], ff['H'][:1], ff['dkx'], ff['dky']) + + +def test_far_to_nearfield_rejects_mismatched_shapes() -> None: + ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8) + bad_h_far = [ff['H'][0], numpy.zeros((8, 4), dtype=complex)] + + with pytest.raises(Exception, match='All fields must be the same shape'): + farfield.far_to_nearfield(ff['E'], bad_h_far, ff['dkx'], ff['dky']) + + +def test_far_to_nearfield_uses_default_and_scalar_padding_shapes() -> None: + ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8) + default_result = farfield.far_to_nearfield( + [field.copy() for field in ff['E']], + [field.copy() for field in ff['H']], + ff['dkx'], + ff['dky'], + ) + scalar_result = farfield.far_to_nearfield( + [field.copy() for field in ff['E']], + [field.copy() for field in ff['H']], + ff['dkx'], + ff['dky'], + padded_size=4, + ) + + assert default_result['E'][0].shape == (8, 8) + assert default_result['H'][0].shape == (8, 8) + assert scalar_result['E'][0].shape == (4, 4) + assert scalar_result['H'][0].shape == (4, 4) diff --git a/meanas/test/test_fdtd_base.py b/meanas/test/test_fdtd_base.py new file mode 100644 index 0000000..41b6ad1 --- /dev/null +++ b/meanas/test/test_fdtd_base.py @@ -0,0 +1,44 @@ +import numpy +from numpy.testing import assert_allclose + +from ..fdmath import functional as fd_functional +from ..fdtd import base + + +DT = 0.25 +SHAPE = (3, 2, 2, 2) +E_FIELD = numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C') / 5.0 +H_FIELD = (numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C') + 1.0) / 7.0 +EPSILON = 1.5 + E_FIELD / 10.0 +MU_FIELD = 2.0 + H_FIELD / 8.0 +MU_SCALAR = 3.0 + + +def test_maxwell_e_without_dxes_matches_unit_spacing_update() -> None: + updater = base.maxwell_e(dt=DT) + expected = E_FIELD + DT * fd_functional.curl_back()(H_FIELD) / EPSILON + + updated = updater(E_FIELD.copy(), H_FIELD.copy(), EPSILON) + + assert_allclose(updated, expected) + + +def test_maxwell_h_without_dxes_and_without_mu_matches_unit_spacing_update() -> None: + updater = base.maxwell_h(dt=DT) + expected = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) + + updated = updater(E_FIELD.copy(), H_FIELD.copy()) + + assert_allclose(updated, expected) + + +def test_maxwell_h_without_dxes_accepts_scalar_and_field_mu() -> None: + updater = base.maxwell_h(dt=DT) + + updated_scalar = updater(E_FIELD.copy(), H_FIELD.copy(), MU_SCALAR) + expected_scalar = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_SCALAR + assert_allclose(updated_scalar, expected_scalar) + + updated_field = updater(E_FIELD.copy(), H_FIELD.copy(), MU_FIELD) + expected_field = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_FIELD + assert_allclose(updated_field, expected_field) diff --git a/meanas/test/test_import_fallbacks.py b/meanas/test/test_import_fallbacks.py new file mode 100644 index 0000000..abee887 --- /dev/null +++ b/meanas/test/test_import_fallbacks.py @@ -0,0 +1,43 @@ +import builtins +import importlib +import pathlib + +import meanas +from ..fdfd import bloch + + +def test_meanas_import_survives_readme_open_failure(monkeypatch) -> None: # type: ignore[no-untyped-def] + original_open = pathlib.Path.open + + def failing_open(self: pathlib.Path, *args, **kwargs): # type: ignore[no-untyped-def] + if self.name == 'README.md': + raise FileNotFoundError('forced README failure') + return original_open(self, *args, **kwargs) + + monkeypatch.setattr(pathlib.Path, 'open', failing_open) + reloaded = importlib.reload(meanas) + + assert reloaded.__version__ == '0.10' + assert reloaded.__author__ == 'Jan Petykiewicz' + assert reloaded.__doc__ is not None + + monkeypatch.undo() + importlib.reload(meanas) + + +def test_bloch_reloads_with_numpy_fft_when_pyfftw_is_unavailable(monkeypatch) -> None: # type: ignore[no-untyped-def] + original_import = builtins.__import__ + + def fake_import(name: str, globals=None, locals=None, fromlist=(), level: int = 0): # type: ignore[no-untyped-def] + if name.startswith('pyfftw'): + raise ImportError('forced pyfftw failure') + return original_import(name, globals, locals, fromlist, level) + + monkeypatch.setattr(builtins, '__import__', fake_import) + reloaded = importlib.reload(bloch) + + assert reloaded.fftn.__module__ == 'numpy.fft' + assert reloaded.ifftn.__module__ == 'numpy.fft' + + monkeypatch.undo() + importlib.reload(bloch) diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index 7bbcd88..d5d3abf 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -162,6 +162,21 @@ def test_waveguide_3d_compute_overlap_e_rejects_empty_overlap_window( ) +def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None: + _epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=2, polarity=1) + + with pytest.raises(ValueError, match='no overlap field support'): + waveguide_3d.compute_overlap_e( + E=numpy.zeros_like(result['E']), + wavenumber=result['wavenumber'], + dxes=dxes, + axis=0, + polarity=1, + slices=slices, + omega=OMEGA, + ) + + def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None: dxes, epsilon, rmin = build_waveguide_cyl_fixture()