diff --git a/docs/user/configuration.md b/docs/user/configuration.md index 7e5f598..21f3944 100644 --- a/docs/user/configuration.md +++ b/docs/user/configuration.md @@ -142,10 +142,16 @@ Common optional keys: - `dark_frame_file` + `blank_frame_file` (both or neither) - `fit.max_points`, `fit.seed`, `fit.loss`, `fit.f_scale`, `fit.max_nfev` -- `fit.initial_C`, `fit.initial_S_i_hat` -- `fit.s3_identifiability_threshold`, `fit.prior_weight`, `fit.c_relative_bounds` +- `fit.initial_C`, `fit.prior_weight`, `fit.c_relative_bounds` - `output_profile`, `output_report`, `output_diagnostics` +Required fit key: + +- `fit.S_i_hat` as an explicit 3-component fixed source state `[S1_hat, S2_hat, S3_hat]` + +`fit.S_i_hat` is treated as a fixed source state during calibration; the +optimiser only fits `C`. Calibration does not provide a default `S_i_hat`. + Load step input notes: - Each `load_steps[i].image_file` may be either: diff --git a/photoelastimetry/calibrate.py b/photoelastimetry/calibrate.py index 96f0f21..115fdfb 100644 --- a/photoelastimetry/calibrate.py +++ b/photoelastimetry/calibrate.py @@ -1,9 +1,9 @@ """ Calibration workflows for photoelastimetry experiments. -This module fits per-wavelength stress-optic coefficients (C), incoming -polarisation state (S_i_hat), and optional detector blank correction from a -multi-load calibration sequence. +This module fits per-wavelength stress-optic coefficients (C) for an explicit +fixed incoming polarisation state (S_i_hat), and optional detector blank +correction from a multi-load calibration sequence. Supported calibration methods: - ``brazilian_disk``: diametrically-loaded disk with analytical stress field. @@ -357,7 +357,7 @@ def validate_calibration_config(config): if n_no_load < 1: warnings.warn( "Calibration usually benefits from at least one no-load step (load ≈ 0); " - "continuing without it will use a default initial S_i_hat guess unless one is provided in fit.initial_S_i_hat.", + "continuing without it may reduce robustness.", stacklevel=2, ) if n_loaded < 3: @@ -368,13 +368,23 @@ def validate_calibration_config(config): ) fit_cfg = dict(config.get("fit", {})) + if "initial_S_i_hat" in fit_cfg: + raise ValueError("fit.initial_S_i_hat is no longer supported; use fit.S_i_hat.") + if "S_i_hat" not in fit_cfg: + raise ValueError( + "Calibration config must include fit.S_i_hat as a 3-component array " "[S1_hat, S2_hat, S3_hat]." + ) + fit_s_i_hat = np.asarray(fit_cfg["S_i_hat"], dtype=float) + if fit_s_i_hat.shape != (3,): + raise ValueError(f"fit.S_i_hat must have shape (3,), got {fit_s_i_hat.shape}.") + fit_cfg.setdefault("max_points", 6000) fit_cfg.setdefault("loss", "soft_l1") fit_cfg.setdefault("f_scale", 0.05) fit_cfg.setdefault("max_nfev", 300) fit_cfg.setdefault("seed", 0) - fit_cfg.setdefault("s3_identifiability_threshold", 0.02) fit_cfg.setdefault("prior_weight", 0.0) + fit_cfg["S_i_hat"] = fit_s_i_hat.tolist() output_profile = config.get("output_profile", "calibration_profile.json5") output_report = config.get("output_report", "calibration_report.md") @@ -739,27 +749,6 @@ def _prepare_sampling_points(roi_mask, max_points, seed): return y_idx, x_idx -def _initial_s_i_hat_from_noload(measured_noload): - """Estimate initial incoming Stokes state from no-load measurements.""" - s1 = float(np.median(measured_noload[..., 0])) - s2 = float(np.median(measured_noload[..., 1])) - magnitude_sq = s1**2 + s2**2 - s3 = float(np.sqrt(max(0.0, 1.0 - magnitude_sq))) - s_i = np.array([s1, s2, s3], dtype=float) - - norm = np.linalg.norm(s_i) - if norm < 1e-12: - return np.array([1.0, 0.0, 0.0], dtype=float) - if norm > 1.0: - s_i = s_i / norm - return s_i - - -def _default_initial_s_i_hat(): - """Return a conservative default incoming Stokes state guess.""" - return np.array([1.0, 0.0, 0.0], dtype=float) - - def _build_dataset(config): """Load calibration images and construct regression dataset.""" steps = config["load_steps"] @@ -814,8 +803,6 @@ def _build_dataset(config): sigma_yy_steps = [] sigma_xy_steps = [] - no_load_tolerance = config["load_zero_tolerance"] - no_load_measurements = [] diagnostic_step_index = int(np.argmax(np.abs(loads))) diagnostic_image = corrected_images[diagnostic_step_index] @@ -845,20 +832,6 @@ def _build_dataset(config): sigma_yy_steps.append(sigma_yy) sigma_xy_steps.append(sigma_xy) - if abs(load) <= no_load_tolerance: - no_load_measurements.append(measured) - - if len(no_load_measurements) == 0: - warnings.warn( - "No no-load measurements were found after processing load steps; " - "using a default initial S_i_hat guess of [1, 0, 0].", - stacklevel=2, - ) - initial_s_i_hat = _default_initial_s_i_hat() - else: - noload = np.concatenate(no_load_measurements, axis=0) - initial_s_i_hat = _initial_s_i_hat_from_noload(noload) - dataset = { "method": config["method"], "wavelengths": config["wavelengths"], @@ -876,7 +849,6 @@ def _build_dataset(config): # Backward-compatible key retained for previous diagnostics consumers. "disk_mask": model_mask, "blank_correction": blank_correction, - "initial_s_i_hat": initial_s_i_hat, "X": X, "Y": Y, "geometry": config["geometry"], @@ -898,29 +870,26 @@ def _normalise_s_i_hat(s_i_hat): return s_i -def _decode_params(params, n_channels, fixed_s3=None): - """Decode parameter vector into C and S_i_hat.""" - C = np.asarray(params[:n_channels], dtype=float) - if fixed_s3 is None: - s_i_hat = np.asarray(params[n_channels : n_channels + 3], dtype=float) - else: - s_i_hat = np.array([params[n_channels], params[n_channels + 1], fixed_s3], dtype=float) - s_i_hat = _normalise_s_i_hat(s_i_hat) - return C, s_i_hat +def _resolve_fixed_s_i_hat(fit_config): + """Resolve the fixed incoming Stokes state used during calibration.""" + s_i_hat = np.asarray(fit_config["S_i_hat"], dtype=float) + if s_i_hat.shape != (3,): + raise ValueError(f"fit.S_i_hat must have shape (3,), got {s_i_hat.shape}.") + return _normalise_s_i_hat(s_i_hat) -def calibration_residuals(params, dataset, fixed_s3=None): +def calibration_residuals(params, dataset, fixed_s_i_hat): """ Compute calibration residuals for least-squares fitting. Parameters ---------- params : ndarray - Parameter vector containing C values and S_i_hat components. + 1D parameter vector containing one C value per wavelength channel. dataset : dict Dataset dictionary from `_build_dataset`. - fixed_s3 : float, optional - If provided, S3 is fixed and only S1/S2 are optimised. + fixed_s_i_hat : array-like + Fixed source state used for the forward model. Returns ------- @@ -928,8 +897,8 @@ def calibration_residuals(params, dataset, fixed_s3=None): 1D residual vector. """ wavelengths = dataset["wavelengths"] - n_channels = wavelengths.size - C, s_i_hat = _decode_params(params, n_channels, fixed_s3=fixed_s3) + C = np.asarray(params[: wavelengths.size], dtype=float) + s_i_hat = _normalise_s_i_hat(fixed_s_i_hat) residual_chunks = [] for measured, sigma_xx, sigma_yy, sigma_xy in zip( @@ -959,7 +928,7 @@ def calibration_residuals(params, dataset, fixed_s3=None): def fit_calibration_parameters(dataset, fit_config): """ - Fit C and S_i_hat from calibration dataset. + Fit C from calibration dataset with fixed S_i_hat. Parameters ---------- @@ -971,7 +940,7 @@ def fit_calibration_parameters(dataset, fit_config): Returns ------- dict - Fit results containing estimated C, S_i_hat, metrics, and optimizer state. + Fit results containing estimated C, fixed S_i_hat, metrics, and optimizer state. """ wavelengths = dataset["wavelengths"] n_channels = wavelengths.size @@ -979,13 +948,7 @@ def fit_calibration_parameters(dataset, fit_config): init_c = _as_float_array( fit_config.get("initial_C", np.full(n_channels, 3e-9)), expected_length=n_channels, name="initial_C" ) - init_s = np.asarray(fit_config.get("initial_S_i_hat", dataset["initial_s_i_hat"]), dtype=float) - if init_s.size == 2: - init_s = np.append(init_s, 0.0) - if init_s.size != 3: - raise ValueError(f"fit.initial_S_i_hat must have length 2 or 3, got {init_s.size}.") - - init_s = _normalise_s_i_hat(init_s) + fixed_s_i_hat = _resolve_fixed_s_i_hat(fit_config) c_relative_bounds = fit_config.get("c_relative_bounds") if c_relative_bounds is not None: @@ -1001,70 +964,26 @@ def fit_calibration_parameters(dataset, fit_config): c_lower = np.full(n_channels, 1e-15) c_upper = np.full(n_channels, 1e-4) - x0_full = np.concatenate([init_c, init_s], axis=0) - lb_full = np.concatenate([c_lower, np.full(3, -1.0)]) - ub_full = np.concatenate([c_upper, np.full(3, 1.0)]) + x0 = init_c.copy() prior_weight = float(fit_config.get("prior_weight", 0.0)) - def residual_full(params): - residual = calibration_residuals(params, dataset, fixed_s3=None) + def residual_c_only(params): + residual = calibration_residuals(params, dataset, fixed_s_i_hat=fixed_s_i_hat) if prior_weight > 0: - prior = prior_weight * (params - x0_full) + prior = prior_weight * (params - x0) return np.concatenate([residual, prior], axis=0) return residual - result_full = least_squares( - residual_full, - x0_full, - bounds=(lb_full, ub_full), + result = least_squares( + residual_c_only, + x0, + bounds=(c_lower, c_upper), loss=fit_config["loss"], f_scale=float(fit_config["f_scale"]), max_nfev=int(fit_config["max_nfev"]), ) - - col_norms = ( - np.linalg.norm(result_full.jac, axis=0) if result_full.jac.size > 0 else np.zeros_like(x0_full) - ) - reference_norm = np.median(col_norms[:-1]) if col_norms.size > 1 else 0.0 - s3_ratio = float(col_norms[-1] / max(reference_norm, 1e-12)) if col_norms.size > 0 else 0.0 - - fallback = False - fallback_reason = None - threshold = float(fit_config.get("s3_identifiability_threshold", 0.02)) - - if (not result_full.success) or (s3_ratio < threshold): - fallback = True - fallback_reason = "s3_identifiability" if s3_ratio < threshold else "full_fit_failed" - - if fallback_reason == "s3_identifiability": - # If S3 is unidentifiable, restart fixed-S3 fit from the configured - # initial guess to avoid inheriting unstable full-model updates. - x0_fix = np.concatenate([init_c, init_s[:2]], axis=0) - else: - c_guess, s_guess = _decode_params(result_full.x, n_channels, fixed_s3=None) - x0_fix = np.concatenate([c_guess, s_guess[:2]], axis=0) - lb_fix = np.concatenate([c_lower, np.full(2, -1.0)]) - ub_fix = np.concatenate([c_upper, np.full(2, 1.0)]) - - def residual_fixed(params): - residual = calibration_residuals(params, dataset, fixed_s3=0.0) - if prior_weight > 0: - prior = prior_weight * (params - x0_fix) - return np.concatenate([residual, prior], axis=0) - return residual - - result = least_squares( - residual_fixed, - x0_fix, - bounds=(lb_fix, ub_fix), - loss=fit_config["loss"], - f_scale=float(fit_config["f_scale"]), - max_nfev=int(fit_config["max_nfev"]), - ) - C_fit, S_fit = _decode_params(result.x, n_channels, fixed_s3=0.0) - else: - result = result_full - C_fit, S_fit = _decode_params(result.x, n_channels, fixed_s3=None) + C_fit = np.asarray(result.x[:n_channels], dtype=float) + S_fit = fixed_s_i_hat residual = result.fun fit_metrics = { @@ -1077,9 +996,6 @@ def residual_fixed(params): "n_residuals": int(residual.size), "n_samples": int(dataset["sample_y"].size), "n_load_steps": int(dataset["loads"].size), - "fallback_used": bool(fallback), - "fallback_reason": fallback_reason, - "s3_identifiability_ratio": float(s3_ratio), } return { diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py index fc17091..76a354f 100644 --- a/tests/test_calibrate.py +++ b/tests/test_calibrate.py @@ -9,12 +9,14 @@ from photoelastimetry.image import simulate_four_step_polarimetry -def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=False): +def _make_synthetic_calibration_case( + tmp_path, noise_scale=0.0, inconsistent_shapes=False, source_state=None, load_scale=1.0 +): height, width = 34, 34 thickness = 0.01 wavelengths = np.array([650e-9, 550e-9, 450e-9], dtype=float) c_true = np.array([2.2e-9, 2.5e-9, 2.8e-9], dtype=float) - s_true = np.array([0.82, 0.18, 0.0], dtype=float) + s_true = np.array([0.82, 0.18, 0.0] if source_state is None else source_state, dtype=float) radius_m = 0.008 pixels_per_meter = 1800.0 @@ -51,10 +53,10 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes= (offset + 1.0 / scale)[np.newaxis, np.newaxis, :, :], (height, width, 3, 4) ).copy() - if noisy: + if noise_scale > 0.0: rng = np.random.default_rng(42) - dark += rng.normal(scale=1e-6, size=dark.shape) - blank += rng.normal(scale=1e-6, size=blank.shape) + dark += rng.normal(scale=0.02 * noise_scale, size=dark.shape) + blank += rng.normal(scale=0.02 * noise_scale, size=blank.shape) dark_file = tmp_path / "dark.npy" blank_file = tmp_path / "blank.npy" @@ -88,15 +90,15 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes= corrected[~disk_mask, :, :] = 1.0 raw = corrected / scale[np.newaxis, np.newaxis, :, :] + offset[np.newaxis, np.newaxis, :, :] - if noisy: - raw += rng.normal(scale=5e-5, size=raw.shape) + if noise_scale > 0.0: + raw += rng.normal(scale=noise_scale, size=raw.shape) if inconsistent_shapes and idx == len(loads) - 1: raw = raw[:-1, :, :, :] image_file = tmp_path / f"load_{idx}.npy" io.save_image(str(image_file), raw.astype(np.float32)) - load_steps.append({"load": float(load), "image_file": str(image_file)}) + load_steps.append({"load": float(load_scale * load), "image_file": str(image_file)}) config = { "method": "brazilian_disk", @@ -119,10 +121,10 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes= "loss": "soft_l1", "f_scale": 0.05, "max_nfev": 250, - "prior_weight": 10000.0 if noisy else 0.0, - "c_relative_bounds": [0.5, 2.0] if noisy else None, + "prior_weight": 10000.0 if noise_scale > 0.0 else 0.0, + "c_relative_bounds": [0.5, 2.0] if noise_scale > 0.0 else None, "initial_C": [2.0e-9, 2.0e-9, 2.0e-9], - "initial_S_i_hat": [0.8, 0.1, 0.0], + "S_i_hat": s_true.tolist(), }, "output_profile": str(tmp_path / "calibration_profile.json5"), "output_report": str(tmp_path / "calibration_report.md"), @@ -201,6 +203,7 @@ def test_build_dataset_accepts_raw_recording_directory_load_steps(tmp_path): {"image_file": str(recording_dir), "load": 2.0}, {"image_file": str(recording_dir), "load": 3.0}, ], + "fit": {"S_i_hat": [1.0, 0.0, 0.0]}, } validated = calibrate.validate_calibration_config(config) @@ -211,13 +214,13 @@ def test_build_dataset_accepts_raw_recording_directory_load_steps(tmp_path): assert dataset["measured_steps"][0].shape[-2:] == (3, 2) -def _make_synthetic_coupon_case(tmp_path): +def _make_synthetic_coupon_case(tmp_path, noise_scale=0.0, source_state=None, load_scale=1.0): height, width = 30, 42 thickness = 0.01 coupon_width_m = 0.012 wavelengths = np.array([650e-9, 550e-9, 450e-9], dtype=float) c_true = np.array([2.4e-9, 2.7e-9, 3.1e-9], dtype=float) - s_true = np.array([0.88, 0.12, 0.0], dtype=float) + s_true = np.array([0.88, 0.12, 0.0] if source_state is None else source_state, dtype=float) # [x0, x1, y0, y1] gauge_roi = np.array([10, 32, 8, 22], dtype=int) @@ -251,6 +254,12 @@ def _make_synthetic_coupon_case(tmp_path): dark_file = tmp_path / "coupon_dark.npy" blank_file = tmp_path / "coupon_blank.npy" + + rng = np.random.default_rng(314) + if noise_scale > 0.0: + dark += rng.normal(scale=0.02 * noise_scale, size=dark.shape) + blank += rng.normal(scale=0.02 * noise_scale, size=blank.shape) + io.save_image(str(dark_file), dark.astype(np.float32)) io.save_image(str(blank_file), blank.astype(np.float32)) @@ -273,10 +282,12 @@ def _make_synthetic_coupon_case(tmp_path): corrected[~roi_mask, :, :] = 1.0 raw = corrected / scale[np.newaxis, np.newaxis, :, :] + offset[np.newaxis, np.newaxis, :, :] + if noise_scale > 0.0: + raw += rng.normal(scale=noise_scale, size=raw.shape) image_file = tmp_path / f"coupon_load_{idx}.npy" io.save_image(str(image_file), raw.astype(np.float32)) - load_steps.append({"load": float(load), "image_file": str(image_file)}) + load_steps.append({"load": float(load_scale * load), "image_file": str(image_file)}) config = { "method": "coupon_test", @@ -300,9 +311,9 @@ def _make_synthetic_coupon_case(tmp_path): "f_scale": 0.05, "max_nfev": 200, "initial_C": [2.2e-9, 2.6e-9, 3.0e-9], - "initial_S_i_hat": [0.85, 0.1, 0.0], - "c_relative_bounds": [0.8, 1.2], - "prior_weight": 8000.0, + "S_i_hat": s_true.tolist(), + "c_relative_bounds": [0.8, 1.2] if noise_scale > 0.0 else None, + "prior_weight": 8000.0 if noise_scale > 0.0 else 0.0, }, "output_profile": str(tmp_path / "coupon_calibration_profile.json5"), "output_report": str(tmp_path / "coupon_calibration_report.md"), @@ -312,27 +323,37 @@ def _make_synthetic_coupon_case(tmp_path): def test_calibration_residuals_near_zero_for_noiseless_data(tmp_path): - config, c_true, s_true = _make_synthetic_calibration_case(tmp_path, noisy=False) + config, c_true, s_true = _make_synthetic_calibration_case(tmp_path) validated = calibrate.validate_calibration_config(config) dataset = calibrate._build_dataset(validated) - params = np.concatenate([c_true, s_true], axis=0) - residual = calibrate.calibration_residuals(params, dataset) + residual = calibrate.calibration_residuals(c_true, dataset, fixed_s_i_hat=s_true) assert residual.ndim == 1 assert np.sqrt(np.mean(residual**2)) < 1e-6 -def test_fit_recovers_c_and_s_i_hat_on_noisy_data(tmp_path): - config, c_true, s_true = _make_synthetic_calibration_case(tmp_path, noisy=True) +@pytest.mark.parametrize( + ("noise_scale", "source_state"), + [ + (0.0, [0.82, 0.18, 0.0]), + (5e-5, [0.82, 0.18, 0.0]), + (0.0, [0.2, 0.1, np.sqrt(0.95)]), + (5e-5, [0.2, 0.1, np.sqrt(0.95)]), + ], +) +def test_fit_recovers_c_with_fixed_s_i_hat_from_synthetic_disk_data(tmp_path, noise_scale, source_state): + config, c_true, s_true = _make_synthetic_calibration_case( + tmp_path, noise_scale=noise_scale, source_state=source_state + ) result = calibrate.run_calibration(config) profile = result["profile"] c_fit = np.asarray(profile["C"], dtype=float) s_fit = np.asarray(profile["S_i_hat"], dtype=float) - assert np.allclose(c_fit, c_true, rtol=0.2, atol=0.0) - assert np.allclose(s_fit[:2], s_true[:2], atol=0.12) + assert np.allclose(c_fit, c_true, rtol=0.25, atol=0.0) + assert np.allclose(s_fit, calibrate._normalise_s_i_hat(s_true), atol=1e-12) assert profile["fit_metrics"]["success"] @@ -368,7 +389,7 @@ def test_blank_correction_coefficients_and_application(): def test_validation_errors_for_bad_configuration(tmp_path): - config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False) + config, _, _ = _make_synthetic_calibration_case(tmp_path) no_noload = dict(config) no_noload["load_steps"] = [dict(step, load=abs(step["load"]) + 1.0) for step in config["load_steps"]] @@ -393,9 +414,22 @@ def test_validation_errors_for_bad_configuration(tmp_path): with pytest.raises(ValueError, match="at least one load_step"): calibrate.validate_calibration_config(empty_steps) + missing_s_i_hat = dict(config) + missing_s_i_hat["fit"] = dict(config["fit"]) + missing_s_i_hat["fit"].pop("S_i_hat") + with pytest.raises(ValueError, match="fit.S_i_hat"): + calibrate.validate_calibration_config(missing_s_i_hat) + + legacy_s_i_hat = dict(config) + legacy_s_i_hat["fit"] = dict(config["fit"]) + legacy_s_i_hat["fit"].pop("S_i_hat") + legacy_s_i_hat["fit"]["initial_S_i_hat"] = [0.8, 0.1, 0.0] + with pytest.raises(ValueError, match="initial_S_i_hat"): + calibrate.validate_calibration_config(legacy_s_i_hat) -def test_build_dataset_without_noload_uses_default_initial_guess(tmp_path): - config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False) + +def test_validation_without_noload_still_accepts_explicit_s_i_hat(tmp_path): + config, _, _ = _make_synthetic_calibration_case(tmp_path) config["load_steps"] = [ dict(step, load=float(index + 1)) for index, step in enumerate(config["load_steps"][:3]) ] @@ -403,14 +437,16 @@ def test_build_dataset_without_noload_uses_default_initial_guess(tmp_path): with pytest.warns(UserWarning, match="no-load"): validated = calibrate.validate_calibration_config(config) - with pytest.warns(UserWarning, match="default initial S_i_hat"): - dataset = calibrate._build_dataset(validated) + dataset = calibrate._build_dataset(validated) - assert np.allclose(dataset["initial_s_i_hat"], np.array([1.0, 0.0, 0.0])) + assert dataset["loads"].shape == (3,) + assert np.allclose( + np.asarray(validated["fit"]["S_i_hat"], dtype=float), np.asarray(config["fit"]["S_i_hat"]) + ) def test_validation_error_for_inconsistent_image_shapes(tmp_path): - config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=True) + config, _, _ = _make_synthetic_calibration_case(tmp_path, inconsistent_shapes=True) validated = calibrate.validate_calibration_config(config) with pytest.raises(ValueError, match="share shape"): @@ -418,7 +454,7 @@ def test_validation_error_for_inconsistent_image_shapes(tmp_path): def test_end_to_end_writes_profile_report_and_diagnostics(tmp_path): - config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False) + config, _, _ = _make_synthetic_calibration_case(tmp_path) result = calibrate.run_calibration(config) assert Path(result["profile_file"]).exists() @@ -459,7 +495,7 @@ def test_coupon_test_end_to_end_fit_and_profile(tmp_path): assert profile["method"] == "coupon_test" assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.2, atol=0.0) - assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float)[:2], s_true[:2], atol=0.12) + assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float), s_true, atol=1e-12) assert profile["fit_metrics"]["success"] loaded = calibrate.load_calibration_profile(result["profile_file"]) @@ -474,3 +510,16 @@ def test_coupon_validation_requires_coupon_geometry_fields(tmp_path): with pytest.raises(ValueError, match="coupon_width_m"): calibrate.validate_calibration_config(bad) + + +@pytest.mark.parametrize("noise_scale", [0.0, 5e-5]) +def test_coupon_calibration_recovers_c_with_fixed_circular_source(tmp_path, noise_scale): + config, c_true, s_true = _make_synthetic_coupon_case( + tmp_path, noise_scale=noise_scale, source_state=[0.0, 0.0, 1.0], load_scale=0.5 + ) + result = calibrate.run_calibration(config) + profile = result["profile"] + + assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.25, atol=0.0) + assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float), s_true, atol=1e-12) + assert profile["fit_metrics"]["success"]