diff --git a/docs/how-it-works/ballistics.md b/docs/how-it-works/ballistics.md index 708f519cb..f57b455fd 100644 --- a/docs/how-it-works/ballistics.md +++ b/docs/how-it-works/ballistics.md @@ -25,14 +25,17 @@ The coefficients depend on the **spin parameter** $S_p = r\omega / v$ — the ratio of surface speed to translational speed: $$ -C_d = C_{d,\text{base}} + k_d S_p +C_d = a + b\,S_p + c\,S_p^2 \qquad -C_l = \frac{C_{l,\text{sat}} \, S_p}{C_{l,\text{half}} + S_p} +C_l = \max\!\left(0,\; d + e\,S_p + f\,S_p^2\right) $$ -Drag rises linearly with spin. Lift follows a Hill-type saturating form: it -approaches a ceiling as spin increases rather than growing without bound, which -is what the wind-tunnel data shows. +Both are second-order polynomials in $S_p$, the form **Ferguson, McNally & +McPhee (2022)** fitted to 1040 measured shots; the coefficients are their +published values. Lift rises with spin, peaks near $S_p \approx 0.52$ and then +falls, which is what the measured data shows and what a saturating form cannot +represent. Above $S_p = 0.75$, the top of the fitted range, both curves are held +at their end value rather than extrapolated. These parametric forms are consistent with **Bearman & Harvey (1976)** and **Kensrud & Smith (2018)** for dimpled balls past the drag crisis @@ -54,10 +57,9 @@ matters across a six-second flight. | `BALL_MASS_KG` | 0.04593 | USGA **maximum**-conforming ball, 45.93 g | | `BALL_RADIUS_M` | 0.02135 | 42.7 mm diameter | | `AIR_DENSITY_STD` | 1.225 kg/m³ | Sea level, 15 °C ISA | -| `CD_BASE` | 0.205 | Drag at zero spin | -| `CD_SPIN_COEFF` | 0.18 | Linear drag rise with $S_p$ | -| `CL_SATURATION` | 0.32 | Lift ceiling | -| `CL_HALF_SP` | 0.15 | $S_p$ at half the lift ceiling | +| `CD_POLY` | (0.1304, 0.9287, -0.8259) | $C_d$ polynomial $(a, b, c)$ | +| `CL_POLY` | (0.0504, 1.2031, -1.1490) | $C_l$ polynomial $(d, e, f)$ | +| `SP_FIT_MAX` | 0.75 | Top of the fitted $S_p$ range; curves held beyond it | | `SPIN_DECAY_RATE` | 0.04 /s | ≈4 %/s | | `GRAVITY` | 9.81 m/s² | | | `DT_SECONDS` | 0.002 | 500 Hz integration | diff --git a/docs/reference/constants.md b/docs/reference/constants.md index 0e2702e2a..9156d52bd 100644 --- a/docs/reference/constants.md +++ b/docs/reference/constants.md @@ -91,10 +91,9 @@ references. | `BALL_MASS_KG` | 0.04593 | | `BALL_RADIUS_M` | 0.02135 | | `AIR_DENSITY_STD` | 1.225 | -| `CD_BASE` | 0.205 | -| `CD_SPIN_COEFF` | 0.18 | -| `CL_SATURATION` | 0.32 | -| `CL_HALF_SP` | 0.15 | +| `CD_POLY` | (0.1304, 0.9287, -0.8259) | +| `CL_POLY` | (0.0504, 1.2031, -1.1490) | +| `SP_FIT_MAX` | 0.75 | | `SPIN_DECAY_RATE` | 0.04 | | `GRAVITY` | 9.81 | | `DT_SECONDS` | 0.002 | diff --git a/scripts/analysis/sweep_ballistic_coeffs.py b/scripts/analysis/sweep_ballistic_coeffs.py index 35514a5e1..d5e9260f4 100644 --- a/scripts/analysis/sweep_ballistic_coeffs.py +++ b/scripts/analysis/sweep_ballistic_coeffs.py @@ -1,11 +1,9 @@ -"""Tune the four aerodynamic coefficients in +"""Tune the six polynomial aerodynamic coefficients in :mod:`openflight.ballistics` against TrackMan-measured carry. -Parameters swept: - CD_BASE — drag coefficient at zero spin - CD_SPIN_COEFF — slope of Cd(Sp) (linear in spin parameter) - CL_SATURATION — Cl asymptote at high Sp - CL_HALF_SP — Sp at which Cl reaches CL_SATURATION/2 +Parameters swept (Cd = a + b*Sp + c*Sp^2, Cl = d + e*Sp + f*Sp^2): + CD_POLY = (a, b, c) + CL_POLY = (d, e, f) Method: scipy.optimize.differential_evolution (global) followed by Nelder-Mead refinement. Loss is overall RMSE on the TrackMan-inputs @@ -77,23 +75,21 @@ # dimpled golf balls in the post-drag-crisis regime; narrow enough that # differential evolution converges in a few minutes. PARAM_BOUNDS: List[Tuple[float, float]] = [ - (0.16, 0.28), # CD_BASE - (0.00, 0.40), # CD_SPIN_COEFF - (0.18, 0.42), # CL_SATURATION - (0.03, 0.35), # CL_HALF_SP + (0.08, 0.30), # Cd a: drag at zero spin + (0.00, 1.80), # Cd b + (-1.80, 0.00), # Cd c + (0.00, 0.15), # Cl d + (0.40, 2.20), # Cl e + (-2.20, -0.30), # Cl f: negative so lift peaks and turns over ] -PARAM_NAMES = ["CD_BASE", "CD_SPIN_COEFF", "CL_SATURATION", "CL_HALF_SP"] -DEFAULT_COEFFS = ( - bl.CD_BASE, - bl.CD_SPIN_COEFF, - bl.CL_SATURATION, - bl.CL_HALF_SP, -) +PARAM_NAMES = ["cd_a", "cd_b", "cd_c", "cl_d", "cl_e", "cl_f"] +Coeffs = Tuple[float, ...] +DEFAULT_COEFFS: Coeffs = (*bl.CD_POLY, *bl.CL_POLY) @dataclass class FitResult: - coeffs: Tuple[float, float, float, float] + coeffs: Coeffs rmse: float preds: List[float] @@ -123,16 +119,16 @@ def _filter_shots(shots: List[TMShot]) -> List[TMShot]: def simulate_with_coeffs( shots: List[TMShot], - coeffs: Tuple[float, float, float, float], + coeffs: Coeffs, ) -> List[float]: """Monkey-patch the ballistics module constants, run simulate() for every shot, restore the originals on exit. - Relies on ``ballistics._cd`` and ``ballistics._cl`` resolving the - constants at call time from the module's global namespace. + Relies on ``ballistics._cd`` and ``ballistics._cl`` resolving + ``CD_POLY``/``CL_POLY`` at call time from the module's global namespace. """ - saved = (bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP) - bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP = coeffs + saved = (bl.CD_POLY, bl.CL_POLY) + bl.CD_POLY, bl.CL_POLY = tuple(coeffs[:3]), tuple(coeffs[3:]) try: results = [] for s in shots: @@ -140,11 +136,11 @@ def simulate_with_coeffs( results.append(traj.carry_yards) return results finally: - bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP = saved + bl.CD_POLY, bl.CL_POLY = saved def make_loss(shots: List[TMShot], measured: np.ndarray): - """Closure that the optimizer can call with a 4-vector.""" + """Closure that the optimizer can call with a 6-vector.""" def _loss(x: np.ndarray) -> float: preds = simulate_with_coeffs(shots, tuple(x)) return float(np.sqrt(np.mean((np.asarray(preds) - measured) ** 2))) @@ -154,7 +150,7 @@ def _loss(x: np.ndarray) -> float: def evaluate( shots: List[TMShot], measured: np.ndarray, - coeffs: Tuple[float, float, float, float], + coeffs: Coeffs, ) -> FitResult: preds = simulate_with_coeffs(shots, coeffs) rmse = float(np.sqrt(np.mean((np.asarray(preds) - measured) ** 2))) @@ -258,7 +254,7 @@ def write_scatter( preds_default: List[float], preds_fit: List[float], out_path: Path, - fit_coeffs: Tuple[float, float, float, float], + fit_coeffs: Coeffs, ) -> None: import matplotlib @@ -318,10 +314,9 @@ def write_scatter( ax.legend(loc="best", fontsize=8) fig.suptitle( - f"Default vs fit (fit: CD_BASE={fit_coeffs[0]:.4f}, " - f"CD_SPIN_COEFF={fit_coeffs[1]:.4f}, " - f"CL_SATURATION={fit_coeffs[2]:.4f}, " - f"CL_HALF_SP={fit_coeffs[3]:.4f})", + "Default vs fit (fit: " + + ", ".join(f"{n}={v:.4f}" for n, v in zip(PARAM_NAMES, fit_coeffs)) + + ")", fontsize=11, ) fig.tight_layout() @@ -464,10 +459,12 @@ def _logged_loss(x): fit_lines = [ "Optimal coefficients (TM-inputs, RMSE objective):", - f" CD_BASE = {best_x[0]:.5f} (default {DEFAULT_COEFFS[0]:.5f})", - f" CD_SPIN_COEFF = {best_x[1]:.5f} (default {DEFAULT_COEFFS[1]:.5f})", - f" CL_SATURATION = {best_x[2]:.5f} (default {DEFAULT_COEFFS[2]:.5f})", - f" CL_HALF_SP = {best_x[3]:.5f} (default {DEFAULT_COEFFS[3]:.5f})", + *( + f" {name:<5} = {fit:>9.5f} (default {default:>9.5f})" + for name, fit, default in zip(PARAM_NAMES, best_x, DEFAULT_COEFFS) + ), + f" CD_POLY = ({best_x[0]:.4f}, {best_x[1]:.4f}, {best_x[2]:.4f})", + f" CL_POLY = ({best_x[3]:.4f}, {best_x[4]:.4f}, {best_x[5]:.4f})", "", f"Baseline RMSE: {baseline.rmse:.3f} yd", f"Fit RMSE: {fit.rmse:.3f} yd", @@ -533,8 +530,7 @@ def _logged_loss(x): loso_lines.append( f"{ho_disp:24s} {len(test_shots):>6d} {sub_fit_rmse:>9.3f} " f"{test_rmse:>10.3f} {test_bias:>+10.3f} " - f"({sub_best[0]:.3f}, {sub_best[1]:.3f}, " - f"{sub_best[2]:.3f}, {sub_best[3]:.3f})" + "(" + ", ".join(f"{v:.3f}" for v in sub_best) + ")" ) fit_lines.extend(loso_lines) elif args.loso: diff --git a/scripts/analysis/validate_ballistics.py b/scripts/analysis/validate_ballistics.py index b705ad01e..c2d2121e4 100644 --- a/scripts/analysis/validate_ballistics.py +++ b/scripts/analysis/validate_ballistics.py @@ -74,6 +74,10 @@ # → 1.1839 kg/m³. We round to 1.184 for the input. TM_FLAT_AIR_DENSITY = 1.184 +# TrackMan reports apex ("Max Height - Height") in feet while its distance +# columns are yards. The model returns apex in yards, so convert to compare. +YD_TO_FT = 3.0 + # Club-name → ClubType map. The comparison CSV uses normalized names like # "driver", "7-iron", "pw"; the raw TM CSV uses "7 Iron", "Driver", "PW". _CLUB_MAP: Dict[str, ClubType] = { @@ -151,6 +155,10 @@ class TMShot: carry_yards: Optional[float] timestamp: str session: str = "" # session label (typically derived from the source filename) + # TrackMan reports apex as "Max Height - Height" in FEET, unlike the + # yard-denominated distance columns. Kept in feet end to end so the + # numbers stay comparable to the raw export. + apex_feet: Optional[float] = None def _default_session_label(path: Path) -> str: @@ -201,6 +209,7 @@ def load_trackman(path: Path, session: Optional[str] = None) -> List[TMShot]: carry_yards=_to_float(row.get("Carry Flat - Length")), timestamp=row.get("Date", "") or "", session=session_label, + apex_feet=_to_float(row.get("Max Height - Height")), )) return shots @@ -284,6 +293,11 @@ class ValidationRow: measured_carry_yards: float model_carry_yards: float delta_yards: float # model - measured + # Apex, in FEET, and only where the source carries it (the paired + # comparison CSV does not). None means "not measured", never zero. + measured_apex_feet: Optional[float] = None + model_apex_feet: Optional[float] = None + delta_apex_feet: Optional[float] = None def _has_required_inputs(*vals) -> bool: @@ -312,6 +326,8 @@ def validate_tm_inputs( spin_source="measured", ) traj = simulate(conditions, air_density=air_density) + # apex_yards is a height in yards; TrackMan's column is feet. + model_apex_feet = traj.apex_yards * YD_TO_FT session_tag = s.session or "default" out.append(ValidationRow( source="tm", @@ -327,6 +343,11 @@ def validate_tm_inputs( measured_carry_yards=s.carry_yards, model_carry_yards=traj.carry_yards, delta_yards=traj.carry_yards - s.carry_yards, + measured_apex_feet=s.apex_feet, + model_apex_feet=model_apex_feet, + delta_apex_feet=( + None if s.apex_feet is None else model_apex_feet - s.apex_feet + ), )) return out diff --git a/src/openflight/ballistics.py b/src/openflight/ballistics.py index e76d53acf..efe5541e1 100644 --- a/src/openflight/ballistics.py +++ b/src/openflight/ballistics.py @@ -38,31 +38,28 @@ BALL_AREA_M2 = math.pi * BALL_RADIUS_M ** 2 AIR_DENSITY_STD = 1.225 # kg/m³ at sea level, 15 °C ISA -# Cd = CD_BASE + CD_SPIN_COEFF * Sp -# Linear rise with spin parameter Sp = r·ω/v. -# Cl = CL_SATURATION * Sp / (CL_HALF_SP + Sp) -# Hill-type saturating form: Cl → CL_SATURATION as Sp → ∞, -# reaches CL_SATURATION/2 at Sp = CL_HALF_SP. -# These are simple parametric forms consistent with Bearman & Harvey (1976) -# and Kensrud & Smith (2018) for dimpled balls past the drag crisis -# (Re ~ 5e4–2e5), which covers the full range of realistic golf shots. +# Cd = CD_POLY[0] + CD_POLY[1]*Sp + CD_POLY[2]*Sp^2 +# Cl = CL_POLY[0] + CL_POLY[1]*Sp + CL_POLY[2]*Sp^2, clamped at >= 0 +# Second-order polynomials in the spin parameter Sp = r*omega/v, the form +# Ferguson, McNally & McPhee (2022, ISEA 14, doi:10.5703/1288284317493) +# fitted to N=1040 shots (GCQuad launch conditions, FlightScope X3 +# carry/apex/offline, Pro V1, Sp 0.02-0.75). The coefficients below are +# their published values, not re-fitted here. +# Lift peaks near Sp ~ 0.52 and turns over, as in Bearman & Harvey (1976) +# and Lyu, Kensrud & Smith (2018). The previous Hill-type saturating form +# was monotonic by construction, so no choice of its two constants could +# represent that peak, and it sat ~4 yd low on apex across every club. +# Beyond Sp = SP_FIT_MAX both curves are held at their end value rather +# than extrapolating a parabola into a region the fit never saw. # -# Fitted with scripts/analysis/sweep_ballistic_coeffs.py against the committed -# TrackMan capture (session_logs/OpenFlight-Test.Normalized.csv, 24 shots, -# differential evolution + Nelder-Mead, rho=1.184 to match TrackMan "Flat"). -# Overall carry RMSE against TrackMan: 24.52 -> 3.97 yd. The previous -# CL_HALF_SP of 0.15 sat well above the low-spin driver regime (driver -# Sp ~ 0.05-0.08), so for drivers the lift curve never left its low-lift -# regime and driver carry ran ~37 yd short. Irons and wedges (Sp ~ 0.21-0.63) -# were already past CL_HALF_SP and ran long instead, which is why the error -# flipped sign by club. Resulting Cl is ~0.13-0.16 for drivers, rising to -# ~0.22-0.23 for irons and wedges; the iron/wedge values sit inside the -# 0.18-0.25 band the cited sources report, while the driver values remain -# below it. -CD_BASE = 0.19071 -CD_SPIN_COEFF = 0.31588 -CL_SATURATION = 0.25544 -CL_HALF_SP = 0.04758 +# Validation, TrackMan measured launch conditions in, simulated flight out, +# Hill form (as re-fit in #229) -> this form: +# committed 24-shot capture: carry RMSE 3.97 -> 2.90 yd +# independent 5-session 279-shot set: carry RMSE 3.48 -> 1.98 yd, +# apex RMSE 4.41 -> 0.97 yd +CD_POLY = (0.1304, 0.9287, -0.8259) +CL_POLY = (0.0504, 1.2031, -1.1490) +SP_FIT_MAX = 0.75 # Exponential spin decay: ω(t) = ω₀·exp(-rate·t). # ~4%/s per Kiratidis & Leinweber (2018); small but matters over ~6 s flights. @@ -166,12 +163,18 @@ def resolve_launch(shot: Shot) -> Optional[LaunchConditions]: ) +def _poly(coeffs: tuple, sp: float) -> float: + return coeffs[0] + coeffs[1] * sp + coeffs[2] * sp * sp + + def _cd(sp: float) -> float: - return CD_BASE + CD_SPIN_COEFF * sp + return _poly(CD_POLY, min(sp, SP_FIT_MAX)) def _cl(sp: float) -> float: - return CL_SATURATION * sp / (CL_HALF_SP + sp) if sp > 0 else 0.0 + if sp <= 0: + return 0.0 + return max(0.0, _poly(CL_POLY, min(sp, SP_FIT_MAX))) def _derivatives( diff --git a/tests/test_ballistics.py b/tests/test_ballistics.py index 635135d18..f1275a4e6 100644 --- a/tests/test_ballistics.py +++ b/tests/test_ballistics.py @@ -6,9 +6,17 @@ import pytest from openflight.ballistics import ( + BALL_RADIUS_M, + CD_POLY, + CL_POLY, CLUB_TYPICAL_SPIN_RPM, + MPH_TO_MPS, + SP_FIT_MAX, SPIN_DECAY_RATE, LaunchConditions, + _cd, + _cl, + _poly, resolve_launch, simulate, ) @@ -178,3 +186,80 @@ def test_zero_launch_angle_does_not_crash(self): def test_total_distance_includes_rollout(self): traj = simulate(_driver()) assert traj.total_yards > traj.carry_yards + + +class TestAeroCoefficientSafeguards: + """Guards on the Cd/Cl quadratics outside their fitted range. + + The Ferguson coefficients are fitted over Sp <= SP_FIT_MAX. Extrapolated + past it both parabolas turn over and cross zero (Cl above Sp ~1.09, Cd + above ~1.25), which would mean lift pulling the ball down and drag + *accelerating* it. Sp is recomputed every integration step from + r*omega/v, and v decays faster than spin, so a normal lob wedge reaches + Sp ~1.5 near apex - inside the negative-drag region. Holding both curves + at their end value is the guard; these tests are what make its removal + fail. + """ + + def test_zero_spin_produces_no_lift(self): + assert _cl(0.0) == 0.0 + + def test_lift_never_negative_above_fit_range(self): + # Raw Cl(1.5) is about -0.73; clamped it must not pull the ball down. + for sp in (SP_FIT_MAX, 1.0, 1.5, 3.0): + assert _cl(sp) >= 0.0, f"Cl({sp}) = {_cl(sp):.4f} is negative" + + def test_drag_stays_positive_above_fit_range(self): + # Raw Cd crosses zero near Sp 1.25; negative drag is unphysical. + for sp in (SP_FIT_MAX, 1.3, 1.5, 3.0): + assert _cd(sp) > 0.0, f"Cd({sp}) = {_cd(sp):.4f} is not positive" + + def test_both_curves_held_flat_above_fit_range(self): + # Beyond the fitted range the value is pinned to the endpoint, so the + # curves are constant there rather than continuing to fall. + for sp in (1.0, 1.5, 3.0): + assert _cd(sp) == pytest.approx(_cd(SP_FIT_MAX)) + assert _cl(sp) == pytest.approx(_cl(SP_FIT_MAX)) + + def test_curves_unclamped_inside_fit_range(self): + # The guard must not disturb the fitted region it sits above. + for sp in (0.1, 0.25, 0.5): + assert _cd(sp) == pytest.approx(_poly(CD_POLY, sp)) + assert _cl(sp) == pytest.approx(_poly(CL_POLY, sp)) + + def test_lob_wedge_reaches_clamped_region_in_flight(self): + """A real lob wedge drives Sp past the fitted range in flight. + + This is why the clamp is not hypothetical: 55 mph / 40 deg / + 10000 rpm peaks near Sp 1.5, where the raw parabola gives + Cd = -0.35 and Cl = -0.73. The clamp is what holds that at the + Sp 0.75 endpoint instead. + + Note this asserts reachability, not the clamp's presence: removing + the clamp *lowers* peak Sp (less drag keeps the ball faster), so it + cannot serve as a red control. The coefficient tests above do that. + """ + cond = LaunchConditions( + ball_speed_mph=55.0, + launch_angle_v=40.0, + launch_angle_h=0.0, + spin_rpm=10000, + spin_axis_deg=0.0, + spin_source="measured", + ) + traj = simulate(cond) + + # Recorded points are downsampled, so this is a lower bound on the + # true peak - which only makes the assertion stricter. + peak_sp = max( + BALL_RADIUS_M * (p.spin_rpm * 2 * math.pi / 60.0) / (p.speed_mph * MPH_TO_MPS) + for p in traj.points + if p.speed_mph > 1e-6 + ) + assert peak_sp > SP_FIT_MAX, ( + f"lob wedge peak Sp {peak_sp:.2f} no longer exceeds {SP_FIT_MAX}; " + f"the clamp would be untested by any realistic shot" + ) + assert 30.0 < traj.carry_yards < 90.0, ( + f"lob wedge carry {traj.carry_yards:.1f} yd is not plausible" + ) diff --git a/tests/test_ballistics_trackman_regression.py b/tests/test_ballistics_trackman_regression.py index 4d301d81c..2aec1b0f6 100644 --- a/tests/test_ballistics_trackman_regression.py +++ b/tests/test_ballistics_trackman_regression.py @@ -36,24 +36,48 @@ TRACKMAN_CSV = _REPO_ROOT / "session_logs" / "OpenFlight-Test.Normalized.csv" -# Per-club RMSE ceilings in yards. Ratcheted down after the aero coefficients -# were re-fit against this same capture (see ballistics.py CD_/CL_ constants). +# Per-club RMSE ceilings in yards. Ratcheted twice: after the aero +# coefficients were re-fit against this capture (#229), and again when the +# Cd/Cl functional form moved to the Ferguson quadratics (#230). # -# before re-fit after re-fit -# club rmse bias rmse bias -# 7-iron 11.64 +11.00 2.96 +1.52 -# driver 38.60 -37.18 1.75 -0.45 -# pitching wedge 13.56 +11.02 6.26 -2.96 -# OVERALL 24.52 -5.05 3.97 -0.45 +# pre-#229 #229 (Hill) quadratic +# club rmse bias rmse bias rmse bias +# 7-iron 11.64 +11.00 2.96 +1.52 1.33 +1.30 +# driver 38.60 -37.18 1.75 -0.45 2.38 -1.33 +# pitching wedge 13.56 +11.02 6.26 -2.96 4.48 +1.83 +# OVERALL 24.52 -5.05 3.97 -0.45 2.90 +0.58 +# +# Driver moved up within its ceiling: the quadratic coefficients are an +# external published fit, not tuned to this capture, and the driver ceiling +# was left where it was rather than raised. # # Ceilings sit slightly above the measured values to absorb float/platform # drift. Lower them again if the model improves; never raise them. RMSE_BUDGET_YARDS = { "driver": 3.0, - "7-iron": 4.0, - "pitching wedge": 7.5, + "7-iron": 2.0, + "pitching wedge": 5.5, +} +OVERALL_RMSE_BUDGET_YARDS = 3.5 + +# Apex ceilings in FEET (TrackMan's "Max Height - Height" column is feet). +# Apex is the strongest evidence the quadratic form is right: nothing was ever +# fitted against it, so the improvement is out-of-sample. +# +# club Hill (#229) quadratic +# 7-iron 18.15 1.53 +# driver 5.08 3.17 +# pitching wedge 18.00 2.81 +# OVERALL 15.05 2.56 +# +# Seeded above measured to absorb float/platform drift, matching the carry +# budgets above. Lower them when the model improves; never raise them. +APEX_RMSE_BUDGET_FEET = { + "driver": 4.0, + "7-iron": 2.0, + "pitching wedge": 3.5, } -OVERALL_RMSE_BUDGET_YARDS = 5.0 +OVERALL_APEX_RMSE_BUDGET_FEET = 3.0 # The reference capture is a fixed, committed file: if the shot count changes, # the fixture changed and every budget above needs re-deriving. @@ -102,3 +126,48 @@ def test_per_club_carry_rmse_within_budget(validation_rows, club): f"(bias {stats['mean']:+.2f}, max |delta| {stats['max_abs']:.2f}, " f"n={stats['n']})." ) + + +def _apex_deltas(rows, club=None): + """Apex errors in feet, skipping rows whose source carried no apex.""" + return [ + r.delta_apex_feet + for r in rows + if r.delta_apex_feet is not None and (club is None or r.club == club) + ] + + +def test_reference_capture_has_apex_measurements(validation_rows): + """Pin apex coverage: a fixture silently losing the column would turn + every apex budget below into a vacuous pass over an empty list.""" + deltas = _apex_deltas(validation_rows) + assert len(deltas) == EXPECTED_SHOT_COUNT, ( + f"{len(deltas)} of {EXPECTED_SHOT_COUNT} reference shots carry an apex " + f"measurement. The fixture's 'Max Height - Height' column changed." + ) + + +def test_overall_apex_rmse_within_budget(validation_rows): + """Model apex vs TrackMan apex, fed TrackMan's own launch conditions.""" + stats = _stats(_apex_deltas(validation_rows)) + assert stats["rmse"] <= OVERALL_APEX_RMSE_BUDGET_FEET, ( + f"Overall apex RMSE {stats['rmse']:.2f} ft exceeds budget " + f"{OVERALL_APEX_RMSE_BUDGET_FEET} ft (bias {stats['mean']:+.2f}, " + f"max |delta| {stats['max_abs']:.2f}, n={stats['n']}). " + f"The ballistic model's flight height regressed." + ) + + +@pytest.mark.parametrize("club", sorted(APEX_RMSE_BUDGET_FEET)) +def test_per_club_apex_rmse_within_budget(validation_rows, club): + """Per-club apex budgets: trajectory shape is club-dependent, and the + quadratic's biggest win (irons and wedges) must not silently erode.""" + deltas = _apex_deltas(validation_rows, club) + assert deltas, f"No reference shots with apex for club {club!r}" + stats = _stats(deltas) + budget = APEX_RMSE_BUDGET_FEET[club] + assert stats["rmse"] <= budget, ( + f"{club}: apex RMSE {stats['rmse']:.2f} ft exceeds budget {budget} ft " + f"(bias {stats['mean']:+.2f}, max |delta| {stats['max_abs']:.2f}, " + f"n={stats['n']})." + )