diff --git a/he11lib/fitting.py b/he11lib/fitting.py index 37bb6f7..e43c3fd 100644 --- a/he11lib/fitting.py +++ b/he11lib/fitting.py @@ -204,23 +204,28 @@ class ModalFitter: def fit_auto( self, planes: list[MeasurementPlane], + camera: CameraModel, + camera_tolerance: CameraModelTolerance, max_order: int = 4, bic_improvement_threshold: float = 10.0, ) -> ReconstructionResult: """Fit with automatic mode-set growth, capped at `max_order`.""" validate_planes(planes) + self._warn_if_degenerate(planes, camera_tolerance) shells = generate_mode_shells(max_order) current_modes = list(shells[0]) - best_result = self.fit(planes, current_modes) - best_bic = self._bic(planes, best_result, current_modes) + best_result = self.fit(planes, current_modes, camera, camera_tolerance) + best_bic = self._bic(planes, best_result, current_modes, camera_tolerance) grew_until_cap = True for shell in shells[1:]: trial_modes = current_modes + shell warm_start = self._warm_start_coefficients(best_result, current_modes) - trial_result = self.fit(planes, trial_modes, initial_coefficients=warm_start) - trial_bic = self._bic(planes, trial_result, trial_modes) + trial_result = self.fit( + planes, trial_modes, camera, camera_tolerance, initial_coefficients=warm_start + ) + trial_bic = self._bic(planes, trial_result, trial_modes, camera_tolerance) if trial_bic < best_bic - bic_improvement_threshold: current_modes = trial_modes @@ -250,10 +255,47 @@ class ModalFitter: coeffs[mode] = amplitude * np.exp(1j * phase) return coeffs - def _bic(self, planes: list[MeasurementPlane], result: ReconstructionResult, modes: list[tuple[int, int]]) -> float: - chi2 = sum(np.sum((r * np.sqrt(self.noise_estimator.weights(p.flux))) ** 2) for r, p in zip(result.residuals, planes)) + def _warn_if_degenerate( + self, planes: list[MeasurementPlane], camera_tolerance: CameraModelTolerance + ) -> None: + """Warn when free camera+z geometry parameters exceed the plane count. + + With only a handful of planes, adding ~7-9 shared camera unknowns + plus one z correction per plane can be practically underdetermined + even though each plane contributes many pixels of data, because + those unknowns are *global* and only weakly constrained by subtle + keystone differences between planes. + """ + free_camera_count = sum(1 for t in tolerance_to_values(camera_tolerance) if t > 0) + free_z_count = sum(1 for p in planes if p.z_tolerance > 0) + free_geometry_count = free_camera_count + free_z_count + + if free_geometry_count > len(planes): + warnings.warn( + f"{free_geometry_count} free camera/z geometry parameters " + f"(from nonzero tolerances) but only {len(planes)} measurement " + "planes; the joint fit may be practically underdetermined. " + "Consider tightening CameraModelTolerance / " + "MeasurementPlane.z_tolerance.", + UserWarning, + stacklevel=3, + ) + + def _bic( + self, + planes: list[MeasurementPlane], + result: ReconstructionResult, + modes: list[tuple[int, int]], + camera_tolerance: CameraModelTolerance, + ) -> float: + chi2 = sum( + np.sum((r * np.sqrt(self.noise_estimator.weights(p.flux))) ** 2) + for r, p in zip(result.residuals, planes) + ) n_data = sum(p.flux.size for p in planes) - n_params = 2 * len(modes) + 4 + free_camera_count = sum(1 for t in tolerance_to_values(camera_tolerance) if t > 0) + free_z_count = sum(1 for p in planes if p.z_tolerance > 0) + n_params = 2 * len(modes) + 4 + free_camera_count + free_z_count return float(chi2 + n_params * np.log(n_data)) def _estimate_uncertainty(self, opt_result, modes, coeffs, total_power): diff --git a/tests/test_fitting.py b/tests/test_fitting.py index cd58cdb..096f32a 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import pytest @@ -244,24 +246,69 @@ def test_fit_recovers_offset_z_within_tolerance(): def test_fit_auto_does_not_add_modes_for_pure_fundamental(): basis = make_basis() - gen = make_generator(basis) + camera = make_camera() + gen = make_generator(basis, camera) planes = gen.generate( - coefficients={(0, 0): 1.0 + 0j}, z_list=Z_LIST, noise_std=1e-4, seed=4 + coefficients={(0, 0): 1.0 + 0j}, z_list=Z_LIST, image_shape=IMAGE_SHAPE, noise_std=1e-4, seed=4 ) fitter = ModalFitter(basis) - result = fitter.fit_auto(planes, max_order=2) + result = fitter.fit_auto(planes, camera=camera, camera_tolerance=zero_tolerance(), max_order=2) assert set(result.purity.keys()) == {(0, 0)} def test_fit_auto_grows_to_include_second_mode(): basis = make_basis() - gen = make_generator(basis) + camera = make_camera() + gen = make_generator(basis, camera) true_coeffs = {(0, 0): 0.9 + 0j, (0, 1): 0.4 + 0j} - planes = gen.generate(coefficients=true_coeffs, z_list=Z_LIST, noise_std=1e-4, seed=5) + planes = gen.generate( + coefficients=true_coeffs, z_list=Z_LIST, image_shape=IMAGE_SHAPE, noise_std=1e-4, seed=5 + ) fitter = ModalFitter(basis) - result = fitter.fit_auto(planes, max_order=2) + result = fitter.fit_auto(planes, camera=camera, camera_tolerance=zero_tolerance(), max_order=2) assert (0, 1) in result.purity or (0, -1) in result.purity + + +def test_fit_auto_warns_when_free_geometry_params_exceed_plane_count(): + basis = make_basis() + camera = make_camera() + gen = make_generator(basis, camera) + planes = gen.generate( + coefficients={(0, 0): 1.0 + 0j}, + z_list=Z_LIST, # 4 planes + image_shape=IMAGE_SHAPE, + z_tolerance=0.05, # +4 free z params + noise_std=1e-4, + seed=11, + ) + + # +7 free camera params (all but the 2 principal_point components) + + # 4 free z params = 11 free geometry params > 4 planes. + generous_tolerance = CameraModelTolerance( + focal_length_px=camera.focal_length_px * 0.05, + position=(0.01, 0.01, 0.01), + orientation_deg=(2.0, 2.0, 2.0), + principal_point=(0.0, 0.0), + ) + + fitter = ModalFitter(basis) + with pytest.warns(UserWarning, match="free camera/z geometry parameters"): + fitter.fit_auto(planes, camera=camera, camera_tolerance=generous_tolerance, max_order=1) + + +def test_fit_auto_does_not_warn_when_geometry_fully_fixed(): + basis = make_basis() + camera = make_camera() + gen = make_generator(basis, camera) + planes = gen.generate( + coefficients={(0, 0): 1.0 + 0j}, z_list=Z_LIST, image_shape=IMAGE_SHAPE, noise_std=1e-4, seed=12 + ) + + fitter = ModalFitter(basis) + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + fitter.fit_auto(planes, camera=camera, camera_tolerance=zero_tolerance(), max_order=1)