From 57d5c7334a4b13d4a0cad00b0b13e9801fdb383f Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Mon, 5 Oct 2026 14:35:44 +0300 Subject: [PATCH] refactor: remove fitter analysis framework Spectrum model, convolution, chi2 fit engine, C(W) and g_V/g_A extractors and result containers are analysis code outside the theory-only scope. Removes the package and its test suites. --- fitter/__init__.py | 38 --- fitter/extractor.py | 373 ---------------------- fitter/fit_engine.py | 541 -------------------------------- fitter/model.py | 238 -------------- fitter/result.py | 132 -------- tests/fitter/test_extractor.py | 161 ---------- tests/fitter/test_fit_engine.py | 176 ----------- tests/fitter/test_model.py | 98 ------ tests/fitter/test_result.py | 100 ------ tests/physics/test_fitter.py | 448 -------------------------- 10 files changed, 2305 deletions(-) delete mode 100644 fitter/__init__.py delete mode 100644 fitter/extractor.py delete mode 100644 fitter/fit_engine.py delete mode 100644 fitter/model.py delete mode 100644 fitter/result.py delete mode 100644 tests/fitter/test_extractor.py delete mode 100644 tests/fitter/test_fit_engine.py delete mode 100644 tests/fitter/test_model.py delete mode 100644 tests/fitter/test_result.py delete mode 100644 tests/physics/test_fitter.py diff --git a/fitter/__init__.py b/fitter/__init__.py deleted file mode 100644 index 4847bf5..0000000 --- a/fitter/__init__.py +++ /dev/null @@ -1,38 +0,0 @@ -""" -Fitter module — main analysis framework for beta spectrometry. - -Provides the analysis workhorse: theoretical model + detector convolution, -multi-parameter fitting, and physics parameter extraction (C(W), g_V, g_A). - -Usage: - from fitter import SpectrumModel, SpectrumFitter, AnalysisResult - - # Build theoretical model with detector response - model = SpectrumModel(spectrum, detector_response, channel_energies) - - # Fit model to experimental data - fitter = SpectrumFitter(model, exp_spectrum) - result = fitter.fit( - x0=[endpoint_keV, normalization, background], - bounds=[(...), (...), (...)], - ) - - # Extract physics parameters - print(result.summary()) -""" - -from .model import SpectrumModel -from .fit_engine import SpectrumFitter, CurveFitter, FitConfig -from .result import AnalysisResult, FitSummary -from .extractor import CWExtractor, GVAExtractor - -__all__ = [ - "SpectrumModel", - "SpectrumFitter", - "CurveFitter", - "FitConfig", - "AnalysisResult", - "FitSummary", - "CWExtractor", - "GVAExtractor", -] diff --git a/fitter/extractor.py b/fitter/extractor.py deleted file mode 100644 index b01925d..0000000 --- a/fitter/extractor.py +++ /dev/null @@ -1,373 +0,0 @@ -""" -C(W) shape factor and g_V/g_A extraction from fitted spectra. - -Provides methods to: -- Extract C(W) from fitted spectrum vs measured spectrum -- Perform Kurie plot analysis for endpoint determination -- Fit C(W) parametrization to extract g_V, g_A - -The beta decay spectrum is: - dΓ/dE = C₀ · p·E·(W₀-E)² · F(Z,E) · C(W) · [1 + δ(E)] - -For allowed transitions with C(W) = 1, the Kurie plot -K(E) = sqrt(dΓ/dE / (p·E·F(Z,E)·[1+δ(E)])) -is linear in E with slope -(W₀-E₀)/E₀. -""" - -from __future__ import annotations - -from dataclasses import dataclass -from typing import TYPE_CHECKING, Any, Optional, Tuple - -import numpy as np - -if TYPE_CHECKING: - from beta_spectrum.spectrum import SpectrumConfig - -from beta_spectrum.utils import T_to_W, momentum - - -@dataclass -class CWExtractionResult: - """Container for C(W) extraction results.""" - - energies_W: np.ndarray - """Total energy grid in m_e units.""" - - energies_keV: np.ndarray - """Kinetic energy grid in keV.""" - - cw_values: np.ndarray - """Extracted C(W) values.""" - - cw_errors: np.ndarray - """Uncertainties on C(W) values.""" - - flux: np.ndarray - """Measured flux (counts per energy bin).""" - - flux_errors: np.ndarray - """Uncertainties on flux.""" - - endpoint_W: float - """Endpoint energy W₀ in m_e units.""" - - endpoint_keV: float - """Endpoint kinetic energy in keV.""" - - fit_result: Optional[Any] = None # FitResult from SpectrumFitter - - def kurie_plot(self) -> Tuple[np.ndarray, np.ndarray]: - """Compute Kurie plot values. - - Returns - ------- - (energies_keV, kurie_values) - """ - p = momentum(self.energies_W) - E = self.energies_W - - # C(W) = (dΓ/dE) / (C₀ · p·E·(W₀-E)² · F · corrections) - # For Kurie plot: K(E) = sqrt(C(W) / (p·E·(W₀-E)²)) - # Simplified: K(E) = sqrt(counts / (p·E·(W₀-E)²)) - phase_space = p * E * (self.endpoint_W - E) ** 2 - mask = phase_space > 0 - - kurie = np.zeros_like(self.energies_keV) - kurie[mask] = np.sqrt( - self.flux[mask] / (phase_space[mask] * self.energies_keV[mask] ** 2) - ) - - return self.energies_keV[mask], kurie[mask] - - -@dataclass -class GVAExtractionResult: - """Container for g_V and g_A extraction results.""" - - gv_eff: float - """Effective vector coupling constant.""" - - gv_error: float - """1σ uncertainty on g_V.""" - - ga_eff: float - """Effective axial coupling constant.""" - - ga_error: float - """1σ uncertainty on g_A.""" - - correlation: float - """Correlation coefficient between g_V and g_A.""" - - chi2_per_dof: float - """Reduced χ².""" - - def summary(self) -> str: - lines = [ - "=" * 60, - "g_V / g_A Extraction", - "=" * 60, - f"g_V^eff = {self.gv_eff:.6f} ± {self.gv_error:.6f}", - f"g_A^eff = {self.ga_eff:.6f} ± {self.ga_error:.6f}", - f"Correlation(g_V, g_A) = {self.correlation:.3f}", - f"χ²/ndof = {self.chi2_per_dof:.3f}", - "=" * 60, - ] - return "\n".join(lines) - - -class CWExtractor: - """Extract C(W) shape factor from fitted spectrum. - - Parameters - ---------- - config : SpectrumConfig - Configuration used for the theoretical model. - """ - - def __init__(self, config: "SpectrumConfig"): - self.config = config - - def extract( - self, - measured_counts: np.ndarray, - measured_errors: np.ndarray, - energies_keV: np.ndarray, - model_counts: np.ndarray, - endpoint_keV: float, - ) -> CWExtractionResult: - """Extract C(W) from measured spectrum and model prediction. - - Parameters - ---------- - measured_counts : np.ndarray - Measured counts per channel. - measured_errors : np.ndarray - Measurement uncertainties. - energies_keV : np.ndarray - Energy grid in keV. - model_counts : np.ndarray - Theoretical model counts (from fit). - endpoint_keV : float - Endpoint energy in keV. - - Returns - ------- - CWExtractionResult - """ - W = T_to_W(energies_keV) - W0 = T_to_W(endpoint_keV) - - # C(W) ∝ measured / model (ratio gives shape factor) - # For allowed transitions, model includes phase space + Fermi + corrections - cw_values = np.where(model_counts > 0, measured_counts / model_counts, 1.0) - cw_errors = np.where( - model_counts > 0, - measured_errors / model_counts, - 0.0, - ) - - # Clip C(W) to reasonable range - cw_values = np.clip(cw_values, 0.01, 10.0) - # Ensure errors are positive - cw_errors = np.maximum(cw_errors, 1e-10) - - return CWExtractionResult( - energies_W=W, - energies_keV=energies_keV, - cw_values=cw_values, - cw_errors=cw_errors, - flux=measured_counts, - flux_errors=measured_errors, - endpoint_W=W0, - endpoint_keV=endpoint_keV, - ) - - def fit_parametrization( - self, - cw_result: "CWExtractionResult", - parametrization: str = "constant", - ) -> np.ndarray: - """Fit C(W) parametrization to extracted values. - - Parameters - ---------- - cw_result : CWExtractionResult - Extracted C(W) values. - parametrization : str - Parametrization type: 'constant', 'linear', 'quadratic'. - - Returns - ------- - np.ndarray - Fitted parameters [a0, a1, ...]. - - Raises - ------ - ValueError - If parametrization is not recognized. - """ - valid = cw_result.energies_W > 0 - W = cw_result.energies_W[valid] - 1.0 - cw = cw_result.cw_values[valid] - - if len(W) < 3: - raise ValueError("Not enough points for fit") - - if parametrization == "constant": - coeffs, cov = np.polyfit(W, cw, 0, cov=True) - return np.array([float(coeffs)]) - elif parametrization == "linear": - coeffs, cov = np.polyfit(W, cw, 1, cov=True) - return coeffs - elif parametrization == "quadratic": - coeffs, cov = np.polyfit(W, cw, 2, cov=True) - return coeffs - else: - raise ValueError(f"Unknown parametrization: {parametrization}") - - def fit_gV_gA( - self, - cw_result: "CWExtractionResult", - M_F: float = 1.0, - M_GT: float = 0.0, - ) -> "GVAExtractionResult": - """Extract effective g_V and g_A from C(W) values. - - Parameters - ---------- - cw_result : CWExtractionResult - Extracted C(W) values. - M_F : float - Fermi matrix element. - M_GT : float - Gamow-Teller matrix element. - - Returns - ------- - GVAExtractionResult - """ - return GVAExtractor.extract_from_curve( - cw_values=cw_result.cw_values, - energies_W=cw_result.energies_W, - endpoint_W=cw_result.endpoint_W, - ) - - def kurie_analysis( - self, - endpoint_keV: float, - ) -> Tuple["CWExtractionResult", Optional[Any]]: - """Perform Kurie plot analysis to extract C(W). - - Placeholder method: in a full implementation, this would - use the fitted model to compute C(W) from experimental data. - - Returns - ------- - (CWExtractionResult, fit_result) - Extracted C(W) and optional fit result. - """ - # For now, return a dummy result — the real implementation - # would use the fitted spectrum model - energies_keV = np.linspace(5, endpoint_keV * 0.98, 50) - W = T_to_W(energies_keV / 1000.0) - cw_values = np.ones_like(W) - cw_errors = np.ones_like(W) * 0.05 - flux = np.ones_like(W) * 10.0 - flux_errors = np.ones_like(W) * 3.0 - W0 = T_to_W(endpoint_keV / 1000.0) - - result = CWExtractionResult( - energies_W=W, - energies_keV=energies_keV, - cw_values=cw_values, - cw_errors=cw_errors, - flux=flux, - flux_errors=flux_errors, - endpoint_W=W0, - endpoint_keV=endpoint_keV, - ) - return result, None - - -class GVAExtractor: - """Extract effective g_V and g_A from C(W) shape factor. - - For allowed transitions: - C(W) ≈ 1 + (g_A/g_V)² · (corrections) - - The ratio g_A^eff/g_V^eff can be extracted from the shape of C(W). - """ - - @staticmethod - def extract_from_curve( - cw_values: np.ndarray, - energies_W: np.ndarray, - endpoint_W: float, - transition_type: str = "F", - ) -> GVAExtractionResult: - """Extract g_V and g_A from C(W) values. - - Parameters - ---------- - cw_values : np.ndarray - Extracted C(W) values. - energies_W : np.ndarray - Total energy grid. - endpoint_W : float - Endpoint energy W₀. - transition_type : str - Transition type (F, GT, etc.). - - Returns - ------- - GVAExtractionResult - """ - # For allowed transitions, C(W) ≈ 1 + b·(W - W₀) where b depends on g_A/g_V - # Simple linear fit of C(W) vs energy near endpoint - mask = (energies_W > 0) & (energies_W < endpoint_W) & (cw_values > 0) - - if np.sum(mask) < 3: - return GVAExtractionResult( - gv_eff=1.0, - gv_error=0.0, - ga_eff=0.0, - ga_error=0.0, - correlation=0.0, - chi2_per_dof=0.0, - ) - - # Fit C(W) = a + b·(W - W₀) - x = energies_W[mask] - endpoint_W - y = cw_values[mask] - - # Linear fit - coeffs, cov = np.polyfit(x, y, 1, cov=True) - a, b = coeffs - - # For allowed transitions: g_V ≈ 1 (CVC), g_A extracted from slope - gv_eff = 1.0 - gv_error = 0.0 - ga_eff = b * endpoint_W # Simplified relation - ga_error = np.sqrt(cov[1, 1]) * endpoint_W - - correlation = ( - cov[0, 1] / (np.sqrt(cov[0, 0]) * np.sqrt(cov[1, 1])) - if cov[0, 0] * cov[1, 1] > 0 - else 0.0 - ) - - # χ² estimate - y_fit = a + b * x - chi2 = float(np.sum(((y - y_fit) / np.maximum(np.abs(y) * 0.01, 1e-10)) ** 2)) - ndof = max(len(x) - 2, 1) - - return GVAExtractionResult( - gv_eff=gv_eff, - gv_error=gv_error, - ga_eff=ga_eff, - ga_error=ga_error, - correlation=float(correlation), - chi2_per_dof=chi2 / ndof if ndof > 0 else chi2, - ) diff --git a/fitter/fit_engine.py b/fitter/fit_engine.py deleted file mode 100644 index a458436..0000000 --- a/fitter/fit_engine.py +++ /dev/null @@ -1,541 +0,0 @@ -""" -Multi-parameter fitting engine for beta spectrum analysis. - -Orchestrates the fitting of theoretical models to experimental data, -including parameter optimization, covariance estimation, and goodness-of-fit. - -Usage: - fitter = SpectrumFitter(model, exp_spectrum, param_names) - result = fitter.fit(x0=[norm, bg], bounds=([0, 0], [100, 1e6])) - print(result.summary()) -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -from typing import Any, Callable, List, Optional, Sequence, Tuple - -import numpy as np -from scipy.optimize import least_squares - - -@dataclass -class FitConfig: - """Configuration for spectrum fitting.""" - - method: str = "trf" - """Optimization method: 'trf', 'dogbox', 'lm'.""" - - max_nfev: int = 1000 - """Maximum function evaluations.""" - - x_scale: str | float = "jac" - """Parameter scaling.""" - - loss: str = "linear" - """Loss function: 'linear', 'soft_l1', 'huber', 'cauchy'.""" - - f_scale: float = 1.0 - """Robust loss scale parameter.""" - - ftol: float = 1e-10 - """Function value tolerance.""" - - xtol: float = 1e-10 - """Parameter change tolerance.""" - - gtol: float = 1e-10 - """Gradient tolerance.""" - - -@dataclass -class FitResult: - """Container for fit results.""" - - parameters: np.ndarray - """Best-fit parameter values.""" - - covariance: np.ndarray - """Parameter covariance matrix.""" - - chi2: float - """Minimum χ² value.""" - - n_points: int - """Number of data points.""" - - n_free: int - """Number of free parameters.""" - - model_values: np.ndarray - """Model values at best-fit parameters.""" - - residuals: np.ndarray - """Data minus model.""" - - success: bool - """Whether optimization converged.""" - - message: str - """Optimization message.""" - - _residual_fn: Optional[Callable] = field(default=None, repr=False) - """Residual function for profile likelihood.""" - - _x_data: Optional[np.ndarray] = field(default=None, repr=False) - """Independent variable data.""" - - _y_data: Optional[np.ndarray] = field(default=None, repr=False) - """Dependent variable data.""" - - _uncertainties: Optional[np.ndarray] = field(default=None, repr=False) - """Measurement uncertainties.""" - - @property - def chi2_per_dof(self) -> float: - """χ² per degree of freedom.""" - return self.chi2 / max(self.n_free, 1) - - @property - def p_value(self) -> float: - """p-value for the fit (survival function of χ²).""" - from scipy.stats import chi2 - - return float(1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free)) - - @property - def parameters_with_errors(self) -> dict[str, tuple[float, float]]: - """Parameter values with 1σ uncertainties.""" - errors = np.sqrt(np.maximum(np.diag(self.covariance), 0.0)) - return { - f"p{i}": (float(self.parameters[i]), float(errors[i])) - for i in range(len(self.parameters)) - } - - @property - def correlation_matrix(self) -> np.ndarray: - """Parameter correlation matrix.""" - diag = np.sqrt(np.maximum(np.diag(self.covariance), 0.0)) - if np.any(diag == 0): - return np.zeros_like(self.covariance) - return self.covariance / np.outer(diag, diag) - - def summary(self, param_names: Optional[List[str]] = None) -> str: - """Human-readable fit summary.""" - lines = [ - "=" * 60, - "Fit Results", - "=" * 60, - f"Converged: {self.success}", - f"Message: {self.message}", - f"χ² = {self.chi2:.4f}", - f"χ²/ndof = {self.chi2_per_dof:.3f} (ndof = {self.n_free})", - f"p-value = {self.p_value:.4f}", - "-" * 60, - "Parameters:", - ] - - errors = np.sqrt(np.maximum(np.diag(self.covariance), 0.0)) - for i in range(len(self.parameters)): - name = param_names[i] if param_names and i < len(param_names) else f"p{i}" - lines.append( - f" {name:20s} = {self.parameters[i]:14.6f} ± {errors[i]:10.6f}" - ) - - lines.append("-" * 60) - lines.append("Correlation matrix (off-diagonal):") - corr = self.correlation_matrix - n = len(self.parameters) - header = " " + "".join(f"{f'p{i}':>12s}" for i in range(n)) - lines.append(header) - for i in range(n): - name = param_names[i] if param_names and i < len(param_names) else f"p{i}" - row = f" {name:8s}" + "".join( - f"{corr[i, j]:12.3f}" for j in range(n) if j != i - ) - lines.append(row) - - lines.append("=" * 60) - return "\n".join(lines) - - def profile_likelihood( - self, - param_index: int, - n_points: int = 50, - param_range: Optional[Tuple[float, float]] = None, - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """Compute profile likelihood for a single parameter.""" - if self._residual_fn is None: - raise RuntimeError("Profile likelihood not available.") - - if param_range is None: - sigma = np.sqrt(abs(self.covariance[param_index, param_index])) - center = self.parameters[param_index] - param_range = (center - 3 * sigma, center + 3 * sigma) - - param_values = np.linspace(param_range[0], param_range[1], n_points) - chi2_values = np.zeros(n_points) - - for i, pval in enumerate(param_values): - params = list(self.parameters) - params[param_index] = pval - res = self._residual_fn(np.array(params)) - chi2_values[i] = np.sum(res**2) - - min_chi2 = np.min(chi2_values) - normalized = chi2_values - min_chi2 - - return param_values, chi2_values, normalized - - def confidence_interval( - self, - param_index: int, - confidence: float = 0.95, - n_points: int = 100, - ) -> Tuple[float, float]: - """Compute confidence interval for a parameter using Δχ² method.""" - from scipy.stats import chi2 - - delta_chi2 = chi2.ppf(confidence, df=1) - param_values, chi2_values, normalized = self.profile_likelihood( - param_index, n_points=n_points - ) - above = normalized >= delta_chi2 - if not np.any(above): - return ( - float(self.parameters[param_index]), - float(self.parameters[param_index]), - ) - crossings = np.where(np.diff(above.astype(int)))[0] - if len(crossings) >= 2: - lower = param_values[crossings[0]] - upper = param_values[crossings[-1] + 1] - elif len(crossings) == 1: - lower = param_values[crossings[0]] - upper = param_values[-1] - else: - lower = param_values[0] - upper = param_values[-1] - return (float(lower), float(upper)) - - -class SpectrumFitter: - """Multi-parameter fitting engine for beta spectrum analysis. - - Wraps a :class:`SpectrumModel` and experimental data to perform - χ² minimization with configurable loss functions and bounds. - - Parameters - ---------- - model : SpectrumModel - Theoretical model (with optional detector convolution). - exp_energies : np.ndarray - Experimental energy grid (keV). - exp_counts : np.ndarray - Experimental counts per channel. - exp_errors : np.ndarray, optional - Experimental uncertainties. If None, assumed to be √counts. - config : FitConfig, optional - Fitting configuration. - """ - - def __init__( - self, - model: Any, - exp_energies: np.ndarray, - exp_counts: np.ndarray, - exp_errors: Optional[np.ndarray] = None, - config: Optional[FitConfig] = None, - ): - self.model = model - self.exp_energies = np.asarray(exp_energies, dtype=np.float64) - self.exp_counts = np.asarray(exp_counts, dtype=np.float64) - - if exp_errors is not None: - self.exp_errors = np.asarray(exp_errors, dtype=np.float64) - else: - self.exp_errors = np.sqrt(np.maximum(self.exp_counts, 1.0)) - - # Filter out zero-uncertainty points - valid = self.exp_errors > 0 - self.exp_energies = self.exp_energies[valid] - self.exp_counts = self.exp_counts[valid] - self.exp_errors = self.exp_errors[valid] - - self.config = config or FitConfig() - self._n_params: Optional[int] = None - self._last_result: Optional[FitResult] = None - - def _residuals(self, params: np.ndarray) -> np.ndarray: - """Compute weighted residuals.""" - model_values = self.model.evaluate(self.exp_energies, params[0], params[1]) - return (self.exp_counts - model_values) / self.exp_errors - - def _chi2(self, params: np.ndarray) -> float: - """Compute χ².""" - model_values = self.model.evaluate(self.exp_energies, params[0], params[1]) - return float(np.sum(((self.exp_counts - model_values) / self.exp_errors) ** 2)) - - def fit( - self, - x0: Optional[Sequence[float]] = None, - bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None, - param_names: Optional[List[str]] = None, - ) -> FitResult: - """Perform χ² minimization. - - Parameters - ---------- - x0 : sequence of float, optional - Initial parameter values. Auto-estimated if None. - bounds : (lower, upper), optional - Parameter bounds. - param_names : list of str, optional - Parameter names for summary output. - - Returns - ------- - FitResult - """ - if x0 is None: - x0 = [1.0, 0.0] - - x0_array = np.asarray(x0, dtype=np.float64) - self._n_params = len(x0_array) - - bounds_array: Optional[Tuple[np.ndarray, np.ndarray]] = None - if bounds is not None: - bounds_array = ( - np.asarray(bounds[0], dtype=np.float64), - np.asarray(bounds[1], dtype=np.float64), - ) - - ls_kwargs: dict[str, Any] = { - "fun": self._residuals, - "x0": x0_array, - "method": self.config.method, - "max_nfev": self.config.max_nfev, - "x_scale": self.config.x_scale, - "loss": self.config.loss, - "f_scale": self.config.f_scale, - "ftol": self.config.ftol, - "xtol": self.config.xtol, - "gtol": self.config.gtol, - } - if bounds is not None: - ls_kwargs["bounds"] = bounds_array - - result = least_squares(**ls_kwargs) - - # Compute covariance - cov = self._compute_covariance(result) - - model_values = self.model.evaluate(self.exp_energies, *result.x) - residuals = self.exp_counts - model_values - - fit_result = FitResult( - parameters=result.x, - covariance=cov, - chi2=self._chi2(result.x), - n_points=len(self.exp_energies), - n_free=self._n_params, - model_values=model_values, - residuals=residuals, - success=result.success, - message=result.message, - _residual_fn=self._residuals, - _x_data=self.exp_energies, - _y_data=self.exp_counts, - _uncertainties=self.exp_errors, - ) - - self._last_result = fit_result - return fit_result - - def _compute_covariance(self, result: Any) -> np.ndarray: - """Estimate parameter covariance from Jacobian.""" - n = len(result.x) - cov = np.zeros((n, n)) - - if result.jac is not None: - jtj = result.jac.T @ result.jac - if np.linalg.cond(jtj) < 1e12: - cov = np.linalg.inv(jtj) - chi2 = self._chi2(result.x) - ndof = max(len(self.exp_energies) - n, 1) - cov = cov * chi2 / ndof - else: - cov = np.linalg.pinv(jtj) - chi2 = self._chi2(result.x) - ndof = max(len(self.exp_energies) - n, 1) - cov = cov * chi2 / ndof - - return cov - - @property - def parameters(self) -> np.ndarray: - """Current best-fit parameters.""" - if self._last_result is None: - raise RuntimeError("No fit result. Call fit() first.") - return self._last_result.parameters - - @property - def covariance(self) -> np.ndarray: - """Current covariance matrix.""" - if self._last_result is None: - raise RuntimeError("No fit result. Call fit() first.") - return self._last_result.covariance - - -class CurveFitter: - """General-purpose χ² curve fitter for arbitrary callable models. - - This is a thin wrapper around scipy's least_squares that provides - the same API as the original beta_spectrum.fitter.CurveFitter. - - Parameters - ---------- - model : callable - Function f(x, *params) that computes model values. - x_data : np.ndarray - Independent variable data. - y_data : np.ndarray - Dependent variable data. - uncertainties : np.ndarray, optional - Measurement uncertainties. If None, assumed to be 1.0. - config : FitConfig, optional - Fitting configuration. - - Examples - -------- - >>> def linear(x, a, b): - ... return a + b * x - >>> fitter = CurveFitter(linear, x, y, uncertainties=unc) - >>> result = fitter.fit([0.0, 0.0]) - """ - - def __init__( - self, - model: Callable, - x_data: np.ndarray, - y_data: np.ndarray, - uncertainties: Optional[np.ndarray] = None, - config: Optional[FitConfig] = None, - ): - self.model = model - self.x_data = np.asarray(x_data, dtype=np.float64) - self.y_data = np.asarray(y_data, dtype=np.float64) - - if uncertainties is not None: - self.uncertainties = np.asarray(uncertainties, dtype=np.float64) - else: - self.uncertainties = np.ones_like(self.y_data) - - # Filter out zero/negative uncertainty points - valid = self.uncertainties > 0 - self.x_data = self.x_data[valid] - self.y_data = self.y_data[valid] - self.uncertainties = self.uncertainties[valid] - - self.config = config or FitConfig() - self._last_result: Optional[FitResult] = None - self._n_params: Optional[int] = None - - def _residuals(self, params: np.ndarray) -> np.ndarray: - """Compute weighted residuals.""" - model_values = self.model(self.x_data, *params) - return (self.y_data - model_values) / self.uncertainties - - def _chi2(self, params: np.ndarray) -> float: - """Compute χ².""" - model_values = self.model(self.x_data, *params) - return float(np.sum(((self.y_data - model_values) / self.uncertainties) ** 2)) - - def fit( - self, - x0: Sequence[float], - bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None, - ) -> FitResult: - """Perform χ² minimization. - - Parameters - ---------- - x0 : sequence of float - Initial parameter values. - bounds : (lower, upper), optional - Parameter bounds. - - Returns - ------- - FitResult - """ - x0_array = np.asarray(x0, dtype=np.float64) - self._n_params = len(x0_array) - - ls_kwargs: dict[str, Any] = { - "fun": self._residuals, - "x0": x0_array, - "method": self.config.method, - "max_nfev": self.config.max_nfev, - "x_scale": self.config.x_scale, - "loss": self.config.loss, - "f_scale": self.config.f_scale, - "ftol": self.config.ftol, - "xtol": self.config.xtol, - "gtol": self.config.gtol, - } - if bounds is not None: - ls_kwargs["bounds"] = ( - np.asarray(bounds[0], dtype=np.float64), - np.asarray(bounds[1], dtype=np.float64), - ) - - result = least_squares(**ls_kwargs) - - # Compute covariance - n = len(result.x) - cov = np.zeros((n, n)) - if result.jac is not None: - jtj = result.jac.T @ result.jac - if np.linalg.cond(jtj) < 1e12: - cov = np.linalg.inv(jtj) - else: - cov = np.linalg.pinv(jtj) - chi2 = self._chi2(result.x) - ndof = max(len(self.x_data) - n, 1) - cov = cov * chi2 / ndof - - model_values = self.model(self.x_data, *result.x) - - fit_result = FitResult( - parameters=result.x, - covariance=cov, - chi2=self._chi2(result.x), - n_points=len(self.x_data), - n_free=self._n_params, - model_values=model_values, - residuals=self.y_data - model_values, - success=result.success, - message=result.message, - _residual_fn=self._residuals, - _x_data=self.x_data, - _y_data=self.y_data, - _uncertainties=self.uncertainties, - ) - - self._last_result = fit_result - return fit_result - - @property - def parameters(self) -> np.ndarray: - """Current best-fit parameters.""" - if self._last_result is None: - raise RuntimeError("No fit result. Call fit() first.") - return self._last_result.parameters - - @property - def covariance(self) -> np.ndarray: - """Current covariance matrix.""" - if self._last_result is None: - raise RuntimeError("No fit result. Call fit() first.") - return self._last_result.covariance diff --git a/fitter/model.py b/fitter/model.py deleted file mode 100644 index 6c8f596..0000000 --- a/fitter/model.py +++ /dev/null @@ -1,238 +0,0 @@ -""" -Theoretical model for fitting — wraps BetaSpectrum + DetectorResponse. - -Provides a callable interface that evaluates the theoretical spectrum, -optionally convolved with detector response, on an arbitrary energy grid. - -Usage: - model = SpectrumModel(beta_spectrum, detector_response, channel_energies) - y_model = model(channel_energies, normalization, background) -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING, Any, Optional - -import numpy as np - -if TYPE_CHECKING: - from beta_spectrum.spectrum import BetaSpectrum - from exp_data.spectrum import ExpSpectrum - - -class SpectrumModel: - """Theoretical beta spectrum model for fitting to experimental data. - - Wraps a :class:`beta_spectrum.BetaSpectrum` and optionally a - :class:`beta_spectrum.DetectorResponse` to produce a callable model - that can be fitted against experimental spectra. - - The model computes: - y(E) = A · [S(E) ⊗ R(E)] + B - - where S(E) is the theoretical spectrum, R(E) is the detector response, - ⊗ denotes convolution, A is a normalization factor, and B is a constant - background level. - - Parameters - ---------- - spectrum : BetaSpectrum - Theoretical spectrum calculator. - detector_response : optional - Detector response object exposing ``channel_energies`` and - ``convolve`` (duck-typed: the DetectorResponse module was removed - from beta_spectrum and is pending a restore-or-delete decision). - If None, no convolution is applied. - channel_energies : np.ndarray, optional - Energy values for detector channels (keV). Used to set up - the response if detector_response is None. - """ - - def __init__( - self, - spectrum: "BetaSpectrum", - detector_response: Optional[Any] = None, - channel_energies: Optional[np.ndarray] = None, - ): - self.spectrum = spectrum - self.detector_response = detector_response - self.channel_energies = channel_energies - - # Cache for pre-computed unconvolved spectrum - self._cached_W: Optional[np.ndarray] = None - self._cached_spectral: Optional[np.ndarray] = None - - def evaluate_unconvolved( - self, - energies_keV: Optional[np.ndarray] = None, - ) -> np.ndarray: - """Evaluate the theoretical spectrum (unconvolved). - - Parameters - ---------- - energies_keV : np.ndarray, optional - Energy grid in keV. Uses the spectrum's internal grid if None. - - Returns - ------- - np.ndarray - Theoretical spectrum values. - """ - if energies_keV is not None: - W = self._keV_to_W(energies_keV) - return self.spectrum(W) - w_grid = getattr(self.spectrum, "_W_grid", None) - if w_grid is None: - return np.array([], dtype=np.float64) - return self.spectrum(w_grid) - - def _keV_to_W(self, energies_keV: np.ndarray) -> np.ndarray: - """Convert kinetic energy in keV to total energy in m_e units.""" - return energies_keV / 510.998950 + 1.0 - - def convolve(self, energies_keV: np.ndarray) -> np.ndarray: - """Convolve theoretical spectrum with detector response. - - Parameters - ---------- - energies_keV : np.ndarray - Energy grid in keV for the detector channels. - - Returns - ------- - np.ndarray - Convolved spectrum on the given energy grid. - - Raises - ------ - RuntimeError - If no detector response is available. - """ - if self.detector_response is None: - raise RuntimeError( - "No detector response configured. " - "Provide a DetectorResponse to SpectrumModel." - ) - - # The detector's channel_energies are in m_e units - W = self.detector_response.channel_energies - # Clip W to valid range (W >= 1.0) before spectrum evaluation - W_safe = np.maximum(W, 1.0 + 1e-12) - spectral = self.spectrum(W_safe) - - # Clip negative/NaN values from spectrum (numerical artifacts) - spectral = np.maximum( - np.nan_to_num(spectral, nan=0.0, posinf=0.0, neginf=0.0), 0.0 - ) - - # Convolve: detector_response.convolve(W, spectrum, normalize) - return self.detector_response.convolve(W, spectral, normalize=True) - - def evaluate( - self, - energies_keV: np.ndarray, - normalization: float = 1.0, - background: float = 0.0, - ) -> np.ndarray: - """Evaluate the full theoretical model (with optional convolution). - - Parameters - ---------- - energies_keV : np.ndarray - Energy grid in keV. - normalization : float - Multiplicative normalization factor. - background : float - Constant background level. - - Returns - ------- - np.ndarray - Model prediction: A · (S ⊗ R) + B - """ - if self.detector_response is not None: - spectral = self.convolve(energies_keV) - else: - spectral = self.spectrum(self._keV_to_W(energies_keV)) - - return normalization * spectral + background - - def residual( - self, - params: np.ndarray, - exp_energies: np.ndarray, - exp_counts: np.ndarray, - exp_errors: np.ndarray, - ) -> np.ndarray: - """Compute weighted residuals for fitting. - - Parameters - ---------- - params : np.ndarray - Fit parameters [normalization, background, ...]. - exp_energies : np.ndarray - Experimental energy grid (keV). - exp_counts : np.ndarray - Experimental counts. - exp_errors : np.ndarray - Experimental uncertainties. - - Returns - ------- - np.ndarray - Weighted residuals: (exp - model) / σ - """ - model_values = self.evaluate(exp_energies, params[0], params[1]) - return (exp_counts - model_values) / exp_errors - - @staticmethod - def default_params( - exp_spectrum: "ExpSpectrum", - ) -> tuple[float, float]: - """Estimate default fit parameters from experimental data. - - Parameters - ---------- - exp_spectrum : ExpSpectrum - Experimental spectrum. - - Returns - ------- - (normalization, background) - Initial parameter estimates. - """ - # Background from low-energy sideband - mask_low = exp_spectrum.energies < np.percentile(exp_spectrum.energies, 5) - background = ( - float(np.mean(exp_spectrum.counts[mask_low])) if np.any(mask_low) else 0.0 - ) - - # Normalization from total counts / integral of unconvolved model - total_counts = float(np.sum(exp_spectrum.counts)) - unconvolved = exp_spectrum.counts - background - norm = total_counts / max( - float(np.trapezoid(unconvolved, exp_spectrum.energies)), 1e-10 - ) - - return (norm, background) - - @staticmethod - def default_bounds( - n_channels: int, - ) -> tuple[list[float], list[float]]: - """Estimate default parameter bounds. - - Parameters - ---------- - n_channels : int - Number of energy channels. - - Returns - ------- - (lower_bounds, upper_bounds) - Parameter bounds for [normalization, background]. - """ - return ( - [0.01, 0.0], - [100.0, 1e6], - ) diff --git a/fitter/result.py b/fitter/result.py deleted file mode 100644 index ae9c220..0000000 --- a/fitter/result.py +++ /dev/null @@ -1,132 +0,0 @@ -""" -Analysis result container for beta spectrum fitting. - -Wraps fit results with physics-specific metadata (nuclide, extracted -C(W) values, g_V, g_A, etc.). -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any, Dict, Optional, Tuple - -import numpy as np - -if TYPE_CHECKING: - from beta_spectrum.spectrum import SpectrumConfig - - -@dataclass -class FitSummary: - """Compact fit summary for reporting.""" - - nuclide: str - """Nuclide identifier (e.g. 'Tc99').""" - - converged: bool - """Whether the fit converged.""" - - chi2_per_dof: float - """Reduced χ².""" - - p_value: float - """Fit p-value.""" - - parameters: Dict[str, Tuple[float, float]] - """Parameter names → (value, error).""" - - def __str__(self) -> str: - lines = [ - f"Fit Summary: {self.nuclide}", - f" Converged: {self.converged}", - f" χ²/ndof: {self.chi2_per_dof:.3f}", - f" p-value: {self.p_value:.4f}", - " Parameters:", - ] - for name, (val, err) in self.parameters.items(): - lines.append(f" {name:20s} = {val:12.6f} ± {err:10.6f}") - return "\n".join(lines) - - -@dataclass -class AnalysisResult: - """Complete analysis result from spectrum fitting. - - Parameters - ---------- - fit_result : FitResult - Low-level fit result from SpectrumFitter. - nuclide : str, optional - Nuclide identifier. - config : SpectrumConfig, optional - Original spectrum configuration. - cw_values : np.ndarray, optional - Extracted C(W) values at each energy bin. - cw_errors : np.ndarray, optional - Uncertainties on C(W). - gv_eff : float, optional - Effective vector coupling constant. - ga_eff : float, optional - Effective axial coupling constant. - endpoint_keV : float, optional - Extracted endpoint energy. - metadata : dict, optional - Additional analysis metadata. - """ - - fit_result: Any # FitResult from SpectrumFitter - nuclide: str = "" - config: Optional["SpectrumConfig"] = None - cw_values: Optional[np.ndarray] = None - cw_errors: Optional[np.ndarray] = None - gv_eff: Optional[float] = None - ga_eff: Optional[float] = None - endpoint_keV: Optional[float] = None - metadata: Dict[str, Any] = field(default_factory=dict) - - @property - def summary(self) -> FitSummary: - """Get a FitSummary from this result.""" - fr = self.fit_result - param_names = self.metadata.get("param_names", None) - params = {} - if hasattr(fr, "parameters_with_errors"): - params = fr.parameters_with_errors - elif hasattr(fr, "parameters"): - errors = np.sqrt(np.maximum(np.diag(fr.covariance), 0.0)) - for i in range(len(fr.parameters)): - name = ( - param_names[i] if param_names and i < len(param_names) else f"p{i}" - ) - params[name] = (float(fr.parameters[i]), float(errors[i])) - - return FitSummary( - nuclide=self.nuclide, - converged=fr.success if hasattr(fr, "success") else False, - chi2_per_dof=fr.chi2_per_dof if hasattr(fr, "chi2_per_dof") else 0.0, - p_value=fr.p_value if hasattr(fr, "p_value") else 0.0, - parameters=params, - ) - - def to_dict(self) -> Dict[str, Any]: - """Serialize to dictionary.""" - return { - "nuclide": self.nuclide, - "converged": ( - self.fit_result.success - if hasattr(self.fit_result, "success") - else False - ), - "chi2_per_dof": ( - self.fit_result.chi2_per_dof - if hasattr(self.fit_result, "chi2_per_dof") - else 0.0 - ), - "p_value": ( - self.fit_result.p_value if hasattr(self.fit_result, "p_value") else 0.0 - ), - "gv_eff": self.gv_eff, - "ga_eff": self.ga_eff, - "endpoint_keV": self.endpoint_keV, - "metadata": self.metadata, - } diff --git a/tests/fitter/test_extractor.py b/tests/fitter/test_extractor.py deleted file mode 100644 index 1987274..0000000 --- a/tests/fitter/test_extractor.py +++ /dev/null @@ -1,161 +0,0 @@ -"""Tests for fitter.extractor.""" - -import numpy as np - -from fitter.extractor import ( - CWExtractor, - GVAExtractor, - CWExtractionResult, - GVAExtractionResult, -) - - -class TestCWExtractor: - """Test CWExtractor.""" - - def test_extract(self): - from beta_spectrum import SpectrumConfig - - config = SpectrumConfig( - Z_parent=43, - Z_daughter=44, - A_number=99, - endpoint_MeV=0.294, - ) - extractor = CWExtractor(config) - - energies_keV = np.linspace(0.01, 0.29, 50) - measured = np.linspace(10.0, 0.1, 50) - errors = np.sqrt(measured) - model = np.linspace(9.5, 0.09, 50) - - result = extractor.extract( - measured_counts=measured, - measured_errors=errors, - energies_keV=energies_keV, - model_counts=model, - endpoint_keV=294.0, - ) - - assert isinstance(result, CWExtractionResult) - assert len(result.cw_values) == 50 - assert len(result.energies_keV) == 50 - assert np.isclose(result.endpoint_keV, 294.0) - assert np.all(np.isfinite(result.cw_values)) - - def test_extract_with_zero_model(self): - from beta_spectrum import SpectrumConfig - - config = SpectrumConfig( - Z_parent=43, - Z_daughter=44, - A_number=99, - endpoint_MeV=0.294, - ) - extractor = CWExtractor(config) - - energies_keV = np.linspace(0.01, 0.29, 50) - measured = np.ones(50) * 5.0 - errors = np.ones(50) * 2.0 - model = np.zeros(50) # All zeros - - result = extractor.extract( - measured_counts=measured, - measured_errors=errors, - energies_keV=energies_keV, - model_counts=model, - endpoint_keV=294.0, - ) - - # C(W) should be 1.0 where model is zero (default) - assert np.all(result.cw_values == 1.0) - - -class TestGVAExtractor: - """Test GVAExtractor.""" - - def test_extract_from_curve_basic(self): - energies_W = np.linspace(1.1, 1.57, 100) - endpoint_W = 1.575 # ~294 keV - cw_values = np.ones_like(energies_W) # C(W) = 1 for allowed - - result = GVAExtractor.extract_from_curve( - cw_values=cw_values, - energies_W=energies_W, - endpoint_W=endpoint_W, - ) - - assert isinstance(result, GVAExtractionResult) - assert np.isclose(result.gv_eff, 1.0) - assert result.chi2_per_dof >= 0 - - def test_extract_from_curve_with_slope(self): - energies_W = np.linspace(1.1, 1.57, 100) - endpoint_W = 1.575 - # C(W) with slight slope - cw_values = 1.0 + 0.1 * (energies_W - endpoint_W) - - result = GVAExtractor.extract_from_curve( - cw_values=cw_values, - energies_W=energies_W, - endpoint_W=endpoint_W, - ) - - assert isinstance(result, GVAExtractionResult) - # g_A should be non-zero due to slope - assert abs(result.ga_eff) > 0 - - def test_extract_from_curve_insufficient_data(self): - energies_W = np.array([1.0, 1.5]) - endpoint_W = 1.575 - cw_values = np.array([1.0, 1.0]) - - result = GVAExtractor.extract_from_curve( - cw_values=cw_values, - energies_W=energies_W, - endpoint_W=endpoint_W, - ) - - assert isinstance(result, GVAExtractionResult) - assert result.chi2_per_dof == 0.0 - - def test_gva_summary(self): - result = GVAExtractionResult( - gv_eff=1.0, - gv_error=0.01, - ga_eff=0.5, - ga_error=0.02, - correlation=-0.3, - chi2_per_dof=1.2, - ) - summary = result.summary() - assert "g_V" in summary - assert "g_A" in summary - assert "1.0" in summary - assert "0.5" in summary - assert "1.2" in summary - - -class TestCWExtractionResult: - """Test CWExtractionResult methods.""" - - def test_kurie_plot(self): - energies_W = np.linspace(1.1, 1.57, 100) - endpoint_W = 1.575 - - result = CWExtractionResult( - energies_W=energies_W, - energies_keV=energies_W * 510.998950 - 510.998950, - cw_values=np.ones_like(energies_W), - cw_errors=np.zeros_like(energies_W), - flux=np.ones_like(energies_W) * 10.0, - flux_errors=np.ones_like(energies_W) * 3.0, - endpoint_W=endpoint_W, - endpoint_keV=294.0, - ) - - E_keV, kurie = result.kurie_plot() - assert len(E_keV) > 0 - assert len(kurie) > 0 - assert np.all(np.isfinite(kurie)) - assert np.all(kurie >= 0) diff --git a/tests/fitter/test_fit_engine.py b/tests/fitter/test_fit_engine.py deleted file mode 100644 index 43f89b2..0000000 --- a/tests/fitter/test_fit_engine.py +++ /dev/null @@ -1,176 +0,0 @@ -"""Tests for fitter.fit_engine.""" - -import numpy as np - -from fitter.fit_engine import SpectrumFitter, FitConfig -from fitter.model import SpectrumModel - - -class TestFitConfig: - """Test FitConfig defaults.""" - - def test_defaults(self) -> None: - config = FitConfig() - assert config.method == "trf" - assert config.max_nfev == 1000 - assert config.loss == "linear" - - -class TestSpectrumFitter: - """Test SpectrumFitter.""" - - def _make_fitter(self) -> "SpectrumFitter": - """Create a simple fitter with synthetic data.""" - - # Create a simple linear model - def linear_model( - energies_keV: np.ndarray, norm: float, bg: float - ) -> np.ndarray: - return norm * (1.0 - energies_keV / 100.0) + bg - - # Synthetic data - energies = np.linspace(0.1, 99.9, 200) - counts = 10.0 * (1.0 - energies / 100.0) + 2.0 - counts += np.random.RandomState(42).normal(0, 0.5, len(energies)) - - # Mock model - class MockModel: - def evaluate( - self, energies_keV: np.ndarray, norm: float, bg: float - ) -> np.ndarray: - return linear_model(energies_keV, norm, bg) - - model = MockModel() - fitter = SpectrumFitter( - model, energies, counts, exp_errors=np.ones_like(counts) * 0.5 - ) - return fitter - - def test_fit_basic(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"]) - - assert result.success - assert result.n_points == 200 - assert result.n_free == 2 - assert result.chi2 > 0 - assert result.chi2_per_dof > 0 - assert len(result.parameters) == 2 - - def test_fit_with_bounds(self) -> None: - fitter = self._make_fitter() - result = fitter.fit( - x0=[10.0, 2.0], - bounds=([1.0, 0.0], [20.0, 10.0]), - ) - assert result.success - - def test_fit_returns_residuals(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - - assert len(result.residuals) == 200 - assert len(result.model_values) == 200 - - def test_parameters_property(self) -> None: - fitter = self._make_fitter() - fitter.fit(x0=[10.0, 2.0]) - params = fitter.parameters - assert len(params) == 2 - - def test_covariance_property(self) -> None: - fitter = self._make_fitter() - fitter.fit(x0=[10.0, 2.0]) - cov = fitter.covariance - assert cov.shape == (2, 2) - - def test_chi2_per_dof(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - assert result.chi2_per_dof == result.chi2 / result.n_free - - def test_p_value(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - assert 0.0 <= result.p_value <= 1.0 - - def test_parameters_with_errors(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - params = result.parameters_with_errors - assert len(params) == 2 - for name, (val, err) in params.items(): - assert isinstance(val, float) - assert isinstance(err, float) - assert err >= 0 - - def test_correlation_matrix(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - corr = result.correlation_matrix - assert corr.shape == (2, 2) - # Diagonal should be 1 - np.testing.assert_allclose(np.diag(corr), 1.0) - - def test_summary(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"]) - summary = result.summary(param_names=["norm", "bg"]) - assert "norm" in summary - assert "bg" in summary - assert "Fit Results" in summary - - def test_profile_likelihood(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - param_vals, chi2_vals, norm_vals = result.profile_likelihood(0, n_points=10) - assert len(param_vals) == 10 - assert len(chi2_vals) == 10 - assert len(norm_vals) == 10 - assert np.min(norm_vals) >= 0 # Normalized χ² should be >= 0 - - def test_confidence_interval(self) -> None: - fitter = self._make_fitter() - result = fitter.fit(x0=[10.0, 2.0]) - ci = result.confidence_interval(0, confidence=0.95) - assert len(ci) == 2 - assert ci[0] <= ci[1] - - -class TestSpectrumFitterWithRealModel: - """Test SpectrumFitter with a real SpectrumModel.""" - - def test_fit_with_real_model(self) -> None: - from beta_spectrum import BetaSpectrum, SpectrumConfig - - config = SpectrumConfig( - Z_parent=43, - Z_daughter=44, - A_number=99, - endpoint_MeV=0.294, - ) - spectrum = BetaSpectrum.from_config(config) - model = SpectrumModel(spectrum) - - # Generate synthetic data from the unconvolved model on a keV grid - energies_keV = np.linspace(1.0, 293.0, 200) - true_norm = 5.0 - true_bg = 0.1 - counts = model.evaluate(energies_keV, true_norm, true_bg) - counts = np.maximum(counts, 0) # Clip negatives - counts += np.random.RandomState(42).normal( - 0, np.sqrt(np.maximum(counts, 1)), len(energies_keV) - ) - counts = np.maximum(counts, 0) - - fitter = SpectrumFitter( - model, - energies_keV, - counts, - exp_errors=np.sqrt(np.maximum(counts, 1)), - ) - result = fitter.fit(x0=[1.0, 0.0], param_names=["norm", "bg"]) - - assert result.success - # Fit should complete without error - assert len(result.parameters) == 2 diff --git a/tests/fitter/test_model.py b/tests/fitter/test_model.py deleted file mode 100644 index 02eb0b8..0000000 --- a/tests/fitter/test_model.py +++ /dev/null @@ -1,98 +0,0 @@ -"""Tests for fitter.model.""" - -import numpy as np -import pytest - -from beta_spectrum import BetaSpectrum, SpectrumConfig -from fitter.model import SpectrumModel - - -class TestSpectrumModel: - """Test SpectrumModel.""" - - def _make_model(self) -> tuple[SpectrumModel, SpectrumConfig]: - """Create a test SpectrumModel (no detector response).""" - config = SpectrumConfig( - Z_parent=43, - Z_daughter=44, - A_number=99, - endpoint_MeV=0.294, - ) - spectrum = BetaSpectrum.from_config(config) - - model = SpectrumModel(spectrum) - return model, config - - def test_evaluate_unconvolved(self) -> None: - model, config = self._make_model() - energies_keV = np.linspace(0.01, 0.29, 50) - values = model.evaluate_unconvolved(energies_keV) - assert len(values) == 50 - assert np.all(np.isfinite(values)) - assert np.all(values >= 0) - - def test_evaluate_no_detector(self) -> None: - model, config = self._make_model() - energies_keV = np.linspace(0.01, 0.29, 50) - values = model.evaluate(energies_keV, normalization=1.0, background=0.0) - assert len(values) == 50 - assert np.all(np.isfinite(values)) - - def test_evaluate_with_normalization(self) -> None: - model, config = self._make_model() - energies_keV = np.linspace(0.01, 0.29, 50) - values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0) - values2 = model.evaluate(energies_keV, normalization=2.0, background=0.0) - np.testing.assert_allclose(values2, values1 * 2.0) - - def test_evaluate_with_background(self) -> None: - model, config = self._make_model() - energies_keV = np.linspace(0.01, 0.29, 50) - values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0) - values2 = model.evaluate(energies_keV, normalization=1.0, background=5.0) - np.testing.assert_allclose(values2, values1 + 5.0) - - def test_convolve_requires_detector(self) -> None: - model, _ = self._make_model() - energies_keV = np.linspace(0.01, 0.29, 50) - with pytest.raises(RuntimeError): - model.convolve(energies_keV) - - def test_default_params(self) -> None: - from exp_data.spectrum import ExpSpectrum - - energies = np.linspace(0, 100, 200) - counts = 50 * np.exp(-0.5 * ((energies - 50) / 5) ** 2) + 10 - exp_spec = ExpSpectrum(energies=energies, counts=counts) - - norm, bg = SpectrumModel.default_params(exp_spec) - assert norm > 0 - assert bg >= 0 - - def test_default_bounds(self) -> None: - lower, upper = SpectrumModel.default_bounds(100) - assert len(lower) == 2 - assert len(upper) == 2 - assert lower[0] < upper[0] - assert lower[1] < upper[1] - - -class TestSpectrumModelKeVToW: - """Test energy unit conversion in SpectrumModel.""" - - def test_keV_to_W(self) -> None: - config = SpectrumConfig( - Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 - ) - model = SpectrumModel(BetaSpectrum.from_config(config)) - # 0 keV → W = 1.0 - W = model._keV_to_W(np.array([0.0])) - assert np.isclose(W[0], 1.0) - - # 511 keV → W = 2.0 - W = model._keV_to_W(np.array([510.998950])) - assert np.isclose(W[0], 2.0) - - # 294 keV → W ≈ 1.575 - W = model._keV_to_W(np.array([294.0])) - assert np.isclose(W[0], 294.0 / 510.998950 + 1.0) diff --git a/tests/fitter/test_result.py b/tests/fitter/test_result.py deleted file mode 100644 index 8bdc711..0000000 --- a/tests/fitter/test_result.py +++ /dev/null @@ -1,100 +0,0 @@ -"""Tests for fitter.result.""" - -import numpy as np - -from fitter.result import AnalysisResult, FitSummary - - -class MockFitResult: - """Mock FitResult for testing AnalysisResult.""" - - def __init__(self): - self.success = True - self.chi2 = 12.5 - self.n_points = 100 - self.n_free = 3 - self.parameters = np.array([5.0, 0.1, 2.0]) - self.covariance = np.array( - [ - [0.01, 0.0, 0.0], - [0.0, 0.001, 0.0], - [0.0, 0.0, 0.04], - ] - ) - self.model_values = np.ones(100) - self.residuals = np.zeros(100) - self.message = "Optimization terminated successfully." - - @property - def chi2_per_dof(self): - return self.chi2 / (self.n_points - self.n_free) - - @property - def p_value(self): - from scipy.stats import chi2 - - return 1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free) - - @property - def parameters_with_errors(self): - errors = np.sqrt(np.diag(self.covariance)) - return { - f"p{i}": (float(self.parameters[i]), float(errors[i])) - for i in range(len(self.parameters)) - } - - -class TestFitSummary: - """Test FitSummary.""" - - def test_str(self): - summary = FitSummary( - nuclide="Tc99", - converged=True, - chi2_per_dof=1.234, - p_value=0.0567, - parameters={"norm": (5.0, 0.1), "bg": (0.1, 0.01)}, - ) - text = str(summary) - assert "Tc99" in text - assert "1.234" in text - assert "5.0" in text - assert "Converged: True" in text - - -class TestAnalysisResult: - """Test AnalysisResult.""" - - def test_summary(self): - mock = MockFitResult() - result = AnalysisResult( - fit_result=mock, - nuclide="Tc99", - metadata={"param_names": ["norm", "bg", "endpoint"]}, - ) - summary = result.summary - assert summary.nuclide == "Tc99" - assert summary.converged is True - assert summary.chi2_per_dof > 0 - - def test_to_dict(self): - mock = MockFitResult() - result = AnalysisResult( - fit_result=mock, - nuclide="Tc99", - gv_eff=0.97, - ga_eff=0.57, - endpoint_keV=294.0, - ) - d = result.to_dict() - assert d["nuclide"] == "Tc99" - assert d["gv_eff"] == 0.97 - assert d["ga_eff"] == 0.57 - assert d["endpoint_keV"] == 294.0 - assert d["converged"] is True - - def test_empty_result(self): - mock = MockFitResult() - result = AnalysisResult(fit_result=mock) - summary = result.summary - assert summary.nuclide == "" diff --git a/tests/physics/test_fitter.py b/tests/physics/test_fitter.py deleted file mode 100644 index e6246e2..0000000 --- a/tests/physics/test_fitter.py +++ /dev/null @@ -1,448 +0,0 @@ -# tests/test_fitter.py — Tests for χ² minimization framework - -""" -Why test this? --------------- -The fitter is the core numerical tool for extracting physical parameters -from beta spectra. It must be correct because all downstream results -(g_V, g_A, C(W) extraction) depend on it. - -We test: - 1. Basic χ² minimization on known functions (quadratic, linear). - 2. Covariance matrix estimation (should match input uncertainties). - 3. Confidence interval computation via profile likelihood. - 4. Goodness-of-fit metrics (χ²/ndf, p-value). - 5. Bounds handling and convergence diagnostics. - 6. Robust loss functions (soft_l1, huber) for outlier resistance. -""" - -import numpy as np - -from fitter.fit_engine import CurveFitter, FitConfig - - -class TestBasicFitting: - """Test χ² minimization on simple known functions.""" - - def test_quadratic_fit(self): - """Fit y = a + b·x + c·x² with known coefficients.""" - np.random.seed(42) - x = np.linspace(0, 10, 50) - true_params = [1.0, 2.0, -0.1] - y = true_params[0] + true_params[1] * x + true_params[2] * x**2 - y += np.random.normal(0, 0.5, len(x)) - - def model(x, a, b, c): - return a + b * x + c * x**2 - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(y) * 0.5) - result = fitter.fit([0.0, 0.0, 0.0]) - - for i, name in enumerate(["a", "b", "c"]): - assert ( - abs(result.parameters[i] - true_params[i]) < 0.5 - ), f"Parameter {name}: expected {true_params[i]}, got {result.parameters[i]}" - - assert result.success, "Fit should converge" - assert result.chi2 > 0, "χ² must be positive" - - def test_linear_fit(self): - """Fit y = a + b·x with known coefficients.""" - np.random.seed(123) - x = np.linspace(0, 5, 30) - true_a, true_b = 3.0, -1.5 - y = true_a + true_b * x - y += np.random.normal(0, 0.2, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(y) * 0.2) - result = fitter.fit([0.0, 0.0]) - - assert abs(result.parameters[0] - true_a) < 0.3 - assert abs(result.parameters[1] - true_b) < 0.3 - assert result.success - - def test_constant_fit(self): - """Fit y = a (constant) with known value.""" - np.random.seed(456) - x = np.linspace(0, 10, 40) - true_a = 5.0 - y = np.full_like(x, true_a) - y += np.random.normal(0, 0.3, len(x)) - - def model(x, a): - return np.full_like(x, a) - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3) - result = fitter.fit([0.0]) - - assert abs(result.parameters[0] - true_a) < 0.2 - assert result.success - - -class TestCovarianceEstimation: - """Test that covariance matrix is correctly estimated.""" - - def test_covariance_matches_uncertainty(self): - """For a constant model, parameter error should match input σ.""" - np.random.seed(789) - x = np.linspace(0, 10, 50) - true_a = 5.0 - sigma = 0.5 - y = np.full_like(x, true_a, dtype=float) - y += np.random.normal(0, sigma, len(x)) - - def model(x, a): - return np.full_like(x, a) - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * sigma) - result = fitter.fit([0.0]) - - error = np.sqrt(result.covariance[0, 0]) - # Expected error: σ / sqrt(N) - expected_error = sigma / np.sqrt(len(x)) - assert ( - abs(error - expected_error) < expected_error * 0.3 - ), f"Covariance error {error:.4f} != expected {expected_error:.4f}" - - def test_covariance_matrix_positive_definite(self): - """Covariance matrix must be positive definite.""" - np.random.seed(321) - x = np.linspace(0, 5, 30) - y = 2.0 + 1.5 * x - 0.3 * x**2 - y += np.random.normal(0, 0.2, len(x)) - - def model(x, a, b, c): - return a + b * x + c * x**2 - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2) - result = fitter.fit([0.0, 0.0, 0.0]) - - # Check eigenvalues are non-negative - eigvals = np.linalg.eigvalsh(result.covariance) - assert np.all(eigvals >= -1e-10), f"Covariance not positive definite: {eigvals}" - - def test_correlation_matrix(self): - """Correlation matrix should have 1.0 on diagonal.""" - np.random.seed(654) - x = np.linspace(0, 10, 50) - y = 1.0 + 2.0 * x + 0.1 * x**2 - y += np.random.normal(0, 0.5, len(x)) - - def model(x, a, b, c): - return a + b * x + c * x**2 - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.5) - result = fitter.fit([0.0, 0.0, 0.0]) - - diag = np.diag(result.correlation_matrix) - assert np.allclose( - diag, 1.0, atol=1e-6 - ), f"Diagonal of correlation matrix: {diag}" - - -class TestGoodnessOfFit: - """Test χ²/ndf and p-value calculations.""" - - def test_chi2_per_dof(self): - """χ²/ndof should be positive and reasonable.""" - np.random.seed(111) - x = np.linspace(0, 10, 50) - y = 2.0 * x - y += np.random.normal(0, 0.3, len(x)) - - def model(x, a): - return a * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3) - result = fitter.fit([0.0]) - - assert result.chi2_per_dof > 0 - assert result.n_points == 50 - assert result.n_free == 1 - assert result.n_free == 1 - - def test_p_value_range(self): - """p-value should be in [0, 1].""" - np.random.seed(222) - x = np.linspace(0, 5, 30) - y = 1.0 + x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.0, 0.0]) - - assert 0.0 <= result.p_value <= 1.0 - - def test_p_value_good_fit(self): - """Good fit should have reasonable p-value.""" - np.random.seed(333) - x = np.linspace(0, 10, 100) - y = 3.0 + 1.5 * x - 0.05 * x**2 - y += np.random.normal(0, 0.2, len(x)) - - def model(x, a, b, c): - return a + b * x + c * x**2 - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2) - result = fitter.fit([0.0, 0.0, 0.0]) - - assert ( - result.p_value > 0.001 - ), f"Good fit should have reasonable p-value, got {result.p_value}" - - -class TestBounds: - """Test parameter bounds handling.""" - - def test_lower_bound(self): - """Test that lower bounds are respected.""" - np.random.seed(444) - x = np.linspace(0, 5, 30) - y = 2.0 + 1.0 * x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([1.5, 0.5], bounds=([1.0, 0.0], [10.0, 5.0])) - - assert ( - result.parameters[0] >= 1.0 - ), f"Parameter a should be >= 1.0, got {result.parameters[0]}" - assert ( - result.parameters[1] >= 0.0 - ), f"Parameter b should be >= 0.0, got {result.parameters[1]}" - - def test_upper_bound(self): - """Test that upper bounds are respected.""" - np.random.seed(555) - x = np.linspace(0, 5, 30) - y = 0.5 + 0.5 * x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.5, 0.5], bounds=([0.0, 0.0], [1.0, 1.0])) - - assert result.parameters[0] <= 1.0 - assert result.parameters[1] <= 1.0 - - -class TestRobustLoss: - """Test robust loss functions for outlier resistance.""" - - def test_soft_l1_loss(self): - """soft_l1 loss should handle outliers better than linear.""" - np.random.seed(666) - x = np.linspace(0, 10, 50) - y = 1.0 + 2.0 * x - y[25] = 50.0 # Add outlier - - def model(x, a, b): - return a + b * x - - # Linear loss — outlier pulls fit away - config_linear = FitConfig(loss="linear") - fitter_linear = CurveFitter( - model, x, y, uncertainties=np.ones_like(x) * 1.0, config=config_linear - ) - result_linear = fitter_linear.fit([0.0, 0.0]) - - # soft_l1 loss — outlier has less influence - config_robust = FitConfig(loss="soft_l1") - fitter_robust = CurveFitter( - model, x, y, uncertainties=np.ones_like(x) * 1.0, config=config_robust - ) - result_robust = fitter_robust.fit([0.0, 0.0]) - - # Robust fit should be closer to true slope - assert ( - abs(result_robust.parameters[1] - 2.0) - < abs(result_linear.parameters[1] - 2.0) * 1.5 - ), "Robust loss should handle outliers better" - - -class TestFitResult: - """Test FitResult dataclass methods.""" - - def test_parameters_with_errors(self): - """parameters_with_errors should return dict of (value, error).""" - np.random.seed(777) - x = np.linspace(0, 5, 30) - y = 2.0 + x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.0, 0.0]) - - errors = result.parameters_with_errors - assert len(errors) == 2 - for key, (val, err) in errors.items(): - assert isinstance(val, float) - assert isinstance(err, float) - assert err >= 0 - - def test_summary(self): - """summary() should return a non-empty string.""" - np.random.seed(888) - x = np.linspace(0, 5, 30) - y = 1.0 + 2.0 * x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.0, 0.0]) - - summary = result.summary(param_names=["a", "b"]) - assert "Fit Results" in summary - assert "a" in summary - assert "b" in summary - assert "χ²" in summary - assert "Converged" in summary - - def test_residuals_property(self): - """residuals should be data - model.""" - np.random.seed(999) - x = np.linspace(0, 5, 30) - y = 2.0 + x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.0, 0.0]) - - expected_residuals = y - result.model_values - np.testing.assert_allclose(result.residuals, expected_residuals) - - -class TestFitConfig: - """Test FitConfig options.""" - - def test_custom_config(self): - """Custom FitConfig should be respected.""" - np.random.seed(101) - x = np.linspace(0, 5, 30) - y = 1.0 + x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - config = FitConfig( - method="trf", - max_nfev=500, - loss="linear", - ftol=1e-10, - ) - fitter = CurveFitter( - model, x, y, uncertainties=np.ones_like(x) * 0.1, config=config - ) - result = fitter.fit([0.0, 0.0]) - - assert result.success - - -class TestConfidenceIntervals: - """Test confidence interval computation.""" - - def test_95_percent_ci(self): - """95% CI should contain the true parameter value.""" - np.random.seed(202) - x = np.linspace(0, 10, 60) - true_a, true_b = 2.0, 1.5 - y = true_a + true_b * x - y += np.random.normal(0, 0.3, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3) - result = fitter.fit([0.0, 0.0]) - - ci_a = result.confidence_interval(0, confidence=0.95) - ci_b = result.confidence_interval(1, confidence=0.95) - - assert ci_a[0] <= true_a <= ci_a[1], f"True a={true_a} not in CI {ci_a}" - assert ci_b[0] <= true_b <= ci_b[1], f"True b={true_b} not in CI {ci_b}" - - def test_ci_width_increases_with_confidence(self): - """Higher confidence level → wider interval.""" - np.random.seed(303) - x = np.linspace(0, 5, 40) - y = 1.0 + 2.0 * x - y += np.random.normal(0, 0.2, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2) - result = fitter.fit([0.0, 0.0]) - - ci_90 = result.confidence_interval(0, confidence=0.90) - ci_95 = result.confidence_interval(0, confidence=0.95) - - width_90 = ci_90[1] - ci_90[0] - width_95 = ci_95[1] - ci_95[0] - - assert ( - width_95 > width_90 - ), f"95% CI width {width_95} should be > 90% width {width_90}" - - -class TestEdgeCases: - """Test edge cases and error handling.""" - - def test_zero_uncertainties_filtered(self): - """Points with zero uncertainty should be filtered out.""" - x = np.array([1.0, 2.0, 3.0, 4.0]) - y = np.array([2.0, 4.0, 6.0, 8.0]) - unc = np.array([1.0, 0.0, 1.0, 0.0]) - - def model(x, a): - return a * x - - fitter = CurveFitter(model, x, y, uncertainties=unc) - assert len(fitter.x_data) == 2, "Zero-uncertainty points should be removed" - assert len(fitter.y_data) == 2 - - def test_negative_uncertainties_filtered(self): - """Points with negative uncertainty should be filtered out.""" - x = np.array([1.0, 2.0, 3.0]) - y = np.array([2.0, 4.0, 6.0]) - unc = np.array([1.0, -1.0, 1.0]) - - def model(x, a): - return a * x - - fitter = CurveFitter(model, x, y, uncertainties=unc) - assert len(fitter.x_data) == 2 - - def test_unbounded_parameters(self): - """Unbounded parameters (±inf) should work.""" - np.random.seed(404) - x = np.linspace(0, 5, 30) - y = 1.0 + 2.0 * x - y += np.random.normal(0, 0.1, len(x)) - - def model(x, a, b): - return a + b * x - - fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1) - result = fitter.fit([0.0, 0.0], bounds=([-np.inf, -np.inf], [np.inf, np.inf])) - assert result.success