Update fit_auto for CameraModel and warn on underdetermined geometry fits
fit_auto now threads camera/camera_tolerance through to fit and _bic (whose parameter count must include any free camera/z unknowns). Emits a UserWarning, not an error, when free camera+z geometry parameters exceed the number of measurement planes -- a new documented degeneracy pitfall. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
+49
-7
@@ -204,23 +204,28 @@ class ModalFitter:
|
|||||||
def fit_auto(
|
def fit_auto(
|
||||||
self,
|
self,
|
||||||
planes: list[MeasurementPlane],
|
planes: list[MeasurementPlane],
|
||||||
|
camera: CameraModel,
|
||||||
|
camera_tolerance: CameraModelTolerance,
|
||||||
max_order: int = 4,
|
max_order: int = 4,
|
||||||
bic_improvement_threshold: float = 10.0,
|
bic_improvement_threshold: float = 10.0,
|
||||||
) -> ReconstructionResult:
|
) -> ReconstructionResult:
|
||||||
"""Fit with automatic mode-set growth, capped at `max_order`."""
|
"""Fit with automatic mode-set growth, capped at `max_order`."""
|
||||||
validate_planes(planes)
|
validate_planes(planes)
|
||||||
|
self._warn_if_degenerate(planes, camera_tolerance)
|
||||||
shells = generate_mode_shells(max_order)
|
shells = generate_mode_shells(max_order)
|
||||||
|
|
||||||
current_modes = list(shells[0])
|
current_modes = list(shells[0])
|
||||||
best_result = self.fit(planes, current_modes)
|
best_result = self.fit(planes, current_modes, camera, camera_tolerance)
|
||||||
best_bic = self._bic(planes, best_result, current_modes)
|
best_bic = self._bic(planes, best_result, current_modes, camera_tolerance)
|
||||||
|
|
||||||
grew_until_cap = True
|
grew_until_cap = True
|
||||||
for shell in shells[1:]:
|
for shell in shells[1:]:
|
||||||
trial_modes = current_modes + shell
|
trial_modes = current_modes + shell
|
||||||
warm_start = self._warm_start_coefficients(best_result, current_modes)
|
warm_start = self._warm_start_coefficients(best_result, current_modes)
|
||||||
trial_result = self.fit(planes, trial_modes, initial_coefficients=warm_start)
|
trial_result = self.fit(
|
||||||
trial_bic = self._bic(planes, trial_result, trial_modes)
|
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:
|
if trial_bic < best_bic - bic_improvement_threshold:
|
||||||
current_modes = trial_modes
|
current_modes = trial_modes
|
||||||
@@ -250,10 +255,47 @@ class ModalFitter:
|
|||||||
coeffs[mode] = amplitude * np.exp(1j * phase)
|
coeffs[mode] = amplitude * np.exp(1j * phase)
|
||||||
return coeffs
|
return coeffs
|
||||||
|
|
||||||
def _bic(self, planes: list[MeasurementPlane], result: ReconstructionResult, modes: list[tuple[int, int]]) -> float:
|
def _warn_if_degenerate(
|
||||||
chi2 = sum(np.sum((r * np.sqrt(self.noise_estimator.weights(p.flux))) ** 2) for r, p in zip(result.residuals, planes))
|
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_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))
|
return float(chi2 + n_params * np.log(n_data))
|
||||||
|
|
||||||
def _estimate_uncertainty(self, opt_result, modes, coeffs, total_power):
|
def _estimate_uncertainty(self, opt_result, modes, coeffs, total_power):
|
||||||
|
|||||||
+53
-6
@@ -1,3 +1,5 @@
|
|||||||
|
import warnings
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
import pytest
|
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():
|
def test_fit_auto_does_not_add_modes_for_pure_fundamental():
|
||||||
basis = make_basis()
|
basis = make_basis()
|
||||||
gen = make_generator(basis)
|
camera = make_camera()
|
||||||
|
gen = make_generator(basis, camera)
|
||||||
planes = gen.generate(
|
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)
|
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)}
|
assert set(result.purity.keys()) == {(0, 0)}
|
||||||
|
|
||||||
|
|
||||||
def test_fit_auto_grows_to_include_second_mode():
|
def test_fit_auto_grows_to_include_second_mode():
|
||||||
basis = make_basis()
|
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}
|
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)
|
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
|
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)
|
||||||
|
|||||||
Reference in New Issue
Block a user