from typing import Literal
import numpy as np
from numpy.typing import NDArray
from perda.core_data_structures import SingleRunData, left_join_data_instances
from perda.units import to_seconds
from scipy.optimize import least_squares
from ..vehicle.irl_vehicle import IrlCar
from .kinematics import BODY_SLIP_ANGLE
from .steering import FRONT_BICYCLE_MODEL_STEER_ANGLE
from .yaw_fit_models import (
BicycleCoeffs,
BicyclePoint,
ChassisParams,
FitComparison,
FitContext,
FitSpec,
Param,
ParamLayout,
Prediction,
WindowArrays,
WindowMetrics,
YawFit,
)
# Speed floor for every 1/V term, guarding the coefficient singularity at V=0.
V_FLOOR_MPS = 1.0
# A fit view: the sub-window, the sign its input/output are mirrored by, and
# the id of the sub-window it came from.
FitView = tuple[WindowArrays, float, int]
def _cg_split(irl_car: IrlCar) -> tuple[float, float]:
"""CG-to-front and CG-to-rear axle distances (m).
SAE convention: ``w_distr_front`` is the *front* weight fraction, so the CG
sits ``wheelbase * (1 - w_distr_front)`` behind the front axle.
"""
wheelbase = irl_car.params.wb
cg_to_front = wheelbase * (1.0 - irl_car.params.w_distr_front)
return cg_to_front, wheelbase - cg_to_front
[docs]
def bicycle_coeffs(
V: NDArray[np.float64],
chassis: ChassisParams,
irl_car: IrlCar,
Ca_f: NDArray[np.float64] | None = None,
Ca_r: NDArray[np.float64] | None = None,
) -> BicycleCoeffs:
"""Linear-bicycle (K, T_z, omega_n, zeta) at each speed. Vectorised over ``V``.
Parameters
----------
V : NDArray[np.float64]
Speed(s) (m/s); must all be strictly positive.
chassis : ChassisParams
The three chassis scalars. ``Izz`` is always taken from here.
irl_car : IrlCar
Source of the mass, wheelbase, and CG split.
Ca_f, Ca_r : NDArray[np.float64], optional
Per-sample cornering stiffnesses, overriding ``chassis``. Used by the
Ca-decay scheduler, where stiffness varies sample-by-sample.
Returns
-------
BicycleCoeffs
Every field shaped like ``V``.
Notes
-----
These are the coefficients of the standard 2-DOF linear bicycle transfer
function from front tire steer (rad) to yaw rate (rad/s)::
H(s; V) = K(V) wn(V)^2 (1 + T_z(V) s) / (s^2 + 2 zeta(V) wn(V) s + wn(V)^2)
All four follow analytically from the three chassis scalars plus the known
geometry, which is why those three scalars are the only thing ``fit_yaw``
optimizes. A free-coefficient LPV polynomial was tried first and proved
unidentifiable on near-limit lapping data: the optimizer slid along the
K/T_z product ridge to nonsense (``T_z ~ 2500 s``, ``K -> 0``) while still
scoring 80-90% VAF. Tying every coefficient back to the bicycle, with hard
prior bounds from the car YAML, removes that ridge.
Inside the prior box ``omega_n^2`` is positive for any sensible understeer
car; the clip below only guards the extreme-oversteer combinations the
optimizer may briefly probe.
"""
V = np.asarray(V, dtype=float)
if np.any(V <= 0):
raise ValueError("V must be strictly positive.")
ca_f = chassis.Ca_f if Ca_f is None else Ca_f
ca_r = chassis.Ca_r if Ca_r is None else Ca_r
m, L = irl_car.params.mass, irl_car.params.wb
lf, lr = _cg_split(irl_car)
Izz = chassis.Izz
omega_n_sq = (ca_f * ca_r * L**2) / (m * Izz * V**2) - (lf * ca_f - lr * ca_r) / Izz
omega_n_sq_eff = np.maximum(omega_n_sq, 1e-9)
omega_n = np.sqrt(omega_n_sq_eff)
two_zeta_wn = (ca_f + ca_r) / (m * V) + (lf**2 * ca_f + lr**2 * ca_r) / (Izz * V)
return BicycleCoeffs(
K=(ca_f * ca_r * L) / (m * Izz * V * omega_n_sq_eff),
T_z=(m * lf * V) / (ca_r * L),
omega_n=omega_n,
zeta=two_zeta_wn / (2.0 * omega_n),
)
[docs]
def bicycle_points(
chassis: ChassisParams, irl_car: IrlCar, speeds: tuple[float, ...]
) -> list[BicyclePoint]:
"""Linear-region bicycle coefficients sampled at each of ``speeds``."""
c = bicycle_coeffs(np.asarray(speeds, dtype=float), chassis, irl_car)
return [
BicyclePoint(
V=float(V),
K=float(c.K[i]),
T_z_ms=float(c.T_z[i]) * 1000.0,
wn_hz=float(c.omega_n[i]) / (2 * np.pi),
zeta=float(c.zeta[i]),
)
for i, V in enumerate(speeds)
]
[docs]
def simulate_lpv_yaw(
t: NDArray[np.float64],
u_rad: NDArray[np.float64],
V: NDArray[np.float64],
chassis: ChassisParams,
irl_car: IrlCar,
Ca_f: NDArray[np.float64] | None = None,
Ca_r: NDArray[np.float64] | None = None,
y_init: float | None = None,
) -> NDArray[np.float64]:
"""Sample-by-sample Tustin biquad with V(t)-scheduled coefficients.
At each sample the four LPV coefficients are evaluated at that sample's
speed (and, when the Ca-decay scheduler is active, that sample's cornering
stiffnesses), then converted to a discrete biquad by Tustin substitution
``s = (2/dt)(z-1)/(z+1)``.
Parameters
----------
t : NDArray[np.float64]
Uniformly-sampled time vector (s).
u_rad : NDArray[np.float64]
Front tire steer input (rad).
V : NDArray[np.float64]
Speed at each sample (m/s); clipped at ``V_FLOOR_MPS``.
chassis : ChassisParams
The three chassis scalars.
irl_car : IrlCar
Source of the mass, wheelbase, and CG split.
Ca_f, Ca_r : NDArray[np.float64], optional
Per-sample cornering stiffnesses; see ``bicycle_coeffs``.
y_init : float, optional
Forced initial output (rad/s). The biquad states are set so ``y[0]``
equals it exactly. Passing the *measured* yaw rate at the sub-window
start is the intended path during fitting: the residual at t=0 is then
zero and the optimizer does not have to fight an initial-condition
mismatch on top of the dynamics. Defaults to the model's own steady
state, ``K(V[0]) * u[0]``.
Returns
-------
NDArray[np.float64]
Simulated yaw rate (rad/s), same length as ``t``.
"""
if not (len(t) == len(u_rad) == len(V)):
raise ValueError("t, u_rad, V must have the same length.")
if len(t) < 2:
return np.zeros_like(u_rad, dtype=float)
dt = float(np.median(np.diff(t)))
if dt <= 0:
raise ValueError("Time vector must be increasing.")
u = np.asarray(u_rad, dtype=float)
V_arr = np.maximum(np.asarray(V, dtype=float), V_FLOOR_MPS)
if Ca_f is not None:
Ca_f = np.broadcast_to(np.asarray(Ca_f, dtype=float), V_arr.shape)
if Ca_r is not None:
Ca_r = np.broadcast_to(np.asarray(Ca_r, dtype=float), V_arr.shape)
c = bicycle_coeffs(V_arr, chassis, irl_car, Ca_f=Ca_f, Ca_r=Ca_r)
K, Tz, wn, zt = c.K, c.T_z, c.omega_n, c.zeta
k = 2.0 / dt
a0 = wn * wn
a1 = 2.0 * zt * wn
b0 = K * a0
b1 = K * a0 * Tz
A0 = k * k + a1 * k + a0
b0_z = (b1 * k + b0) / A0
b1_z = (2.0 * b0) / A0
b2_z = (b0 - b1 * k) / A0
a1_z = (-2.0 * k * k + 2.0 * a0) / A0
a2_z = (k * k - a1 * k + a0) / A0
y = np.empty(len(u), dtype=float)
y0 = float(y_init) if y_init is not None else float(K[0] * u[0])
s1 = y0 - b0_z[0] * u[0]
s2 = b2_z[0] * u[0] - a2_z[0] * y0
for i in range(len(u)):
yi = b0_z[i] * u[i] + s1
s1 = b1_z[i] * u[i] - a1_z[i] * yi + s2
s2 = b2_z[i] * u[i] - a2_z[i] * yi
y[i] = yi
return y
[docs]
def apply_relaxation_lag(
u_rad: NDArray[np.float64],
V: NDArray[np.float64],
t: NDArray[np.float64],
L_f: float,
) -> NDArray[np.float64]:
"""Lag the steering input by a first-order LPV filter, ``tau(V) = L_f / V``.
Approximates front-tire relaxation by lagging the input rather than the
full kinematic slip angle: exact at low frequency, slightly off at high
frequency because the body terms (beta, lf*r/V) are not lagged. A cheap
pre-filter that captures relaxation length's dominant phase effect without
a state-space rewrite. ``L_f <= 0`` is a no-op.
Parameters
----------
u_rad : NDArray[np.float64]
Steering input (rad).
V : NDArray[np.float64]
Speed at each sample (m/s); clipped at ``V_FLOOR_MPS``.
t : NDArray[np.float64]
Uniformly-sampled time vector (s).
L_f : float
Front relaxation length (m).
Returns
-------
NDArray[np.float64]
Lagged input, same length as ``u_rad``.
"""
if L_f <= 0:
return u_rad
if not (len(u_rad) == len(V) == len(t)):
raise ValueError("u_rad, V, t must have the same length.")
if len(u_rad) < 2:
return np.asarray(u_rad, dtype=float).copy()
dt = float(np.median(np.diff(t)))
if dt <= 0:
raise ValueError("Time vector must be increasing.")
tau = L_f / np.maximum(np.asarray(V, dtype=float), V_FLOOR_MPS)
a = dt / (2.0 * tau + dt)
u = np.asarray(u_rad, dtype=float)
y = np.empty_like(u)
y[0] = u[0]
for i in range(1, len(u)):
y[i] = (1.0 - 2.0 * a[i]) * y[i - 1] + a[i] * (u[i] + u[i - 1])
return y
[docs]
def axle_slip_angles(
w: WindowArrays, irl_car: IrlCar
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
"""Per-sample front and rear slip-angle magnitudes (rad).
``|alpha_f| = |delta - (beta + lf r / V)|`` and
``|alpha_r| = |-beta + lr r / V|``. Both derive purely from *measured*
signals, so the optimizer sees them as a fixed time-series -- there is no
feedback through the simulated state during fitting.
Parameters
----------
w : WindowArrays
Must carry both ``body_slip_rad`` and ``bicycle_steer_rad``.
irl_car : IrlCar
Source of the wheelbase and CG split.
Returns
-------
tuple[NDArray[np.float64], NDArray[np.float64]]
``(|alpha_f|, |alpha_r|)``, each the length of ``w.t``.
Raises
------
KeyError
If the body-slip or bicycle-steer channel is absent from ``w``.
"""
if w.body_slip_rad is None:
raise KeyError(
f"Ca decay requires the body-slip channel '{BODY_SLIP_ANGLE}'. "
"Re-run segmentation with add_body_slip_angle()."
)
if w.bicycle_steer_rad is None:
raise KeyError(
f"Ca decay requires the bicycle-steer channel "
f"'{FRONT_BICYCLE_MODEL_STEER_ANGLE}'. Re-run segmentation with "
"add_bicycle_steer_angle()."
)
lf, lr = _cg_split(irl_car)
r = w.yaw_rate
V = np.maximum(w.speed, V_FLOOR_MPS)
beta = w.body_slip_rad
delta = w.bicycle_steer_rad
alpha_f = np.abs(delta - (beta + lf * r / V))
alpha_r = np.abs(-beta + lr * r / V)
return alpha_f, alpha_r
[docs]
def build_param_layout(
spec: FitSpec,
irl_car: IrlCar,
bounds_pct: dict[str, float] | None = None,
u_off_max_deg: float = 10.0,
alpha_init: float = 1.0,
alpha_bounds: tuple[float, float] = (0.0, 50.0),
p_init: float = 1.0,
p_bounds: tuple[float, float] = (0.5, 4.0),
L_f_init: float = 0.3,
L_f_bounds: tuple[float, float] = (0.05, 1.5),
pin: dict[str, float] | None = None,
) -> ParamLayout:
"""Build the parameter vector for ``spec``, boxed by the car's priors.
``Ca_f``, ``Ca_r``, ``Izz`` and ``u_off`` are always present. A Ca
scheduler adds ``alpha_f``, ``alpha_r``, ``p``; front relaxation adds
``L_f``. The three chassis priors and their bounds come from the car
YAML's ``dynamics_setup``.
Parameters
----------
spec : FitSpec
Decides which optional parameters are appended.
irl_car : IrlCar
Source of the ``dynamics_setup`` priors.
bounds_pct : dict[str, float], optional
Overrides the YAML's ``prior_bounds``. Keys ``Ca_front_pct``,
``Ca_rear_pct``, ``Izz_pct``.
u_off_max_deg : float, optional
Symmetric bound on the steering-offset parameter (deg).
alpha_init, alpha_bounds : optional
Initial guess and bounds for the per-axle Ca-decay rates.
p_init, p_bounds : optional
Initial guess and bounds for the shared Ca-decay exponent.
L_f_init, L_f_bounds : optional
Initial guess and bounds for the front relaxation length (m).
pin : dict[str, float], optional
Parameters to hold fixed at the given values instead of fitting.
Returns
-------
ParamLayout
Raises
------
ValueError
If ``irl_car`` carries no ``dynamics_setup`` priors, or ``pin`` names a
parameter this spec does not have.
"""
dyn = irl_car.dynamics_setup
if dyn is None:
raise ValueError(
"build_param_layout: this IrlCar has no yaw-response priors. Add a "
"`dynamics_setup` block (Ca_front_N_per_rad, Ca_rear_N_per_rad, "
"Izz_kg_m2, prior_bounds) to the car YAML."
)
b = dyn.prior_bounds
pct = bounds_pct or {
"Ca_front_pct": b.Ca_front_pct,
"Ca_rear_pct": b.Ca_rear_pct,
"Izz_pct": b.Izz_pct,
}
def prior(name: str, value: float, pct_key: str) -> Param:
"""A chassis scalar, boxed at +/- its prior percentage."""
frac = pct[pct_key] / 100.0
return Param(
name=name,
init=value,
lo=value * (1.0 - frac),
hi=value * (1.0 + frac),
scale=abs(value),
)
params = [
prior("Ca_f", dyn.Ca_front_N_per_rad, "Ca_front_pct"),
prior("Ca_r", dyn.Ca_rear_N_per_rad, "Ca_rear_pct"),
prior("Izz", dyn.Izz_kg_m2, "Izz_pct"),
]
if spec.ca_scheduler != "none":
lo, hi = alpha_bounds
params += [
Param(name="alpha_f", init=alpha_init, lo=lo, hi=hi, scale=1.0),
Param(name="alpha_r", init=alpha_init, lo=lo, hi=hi, scale=1.0),
Param(name="p", init=p_init, lo=p_bounds[0], hi=p_bounds[1], scale=1.0),
]
if spec.relaxation == "front":
params.append(
Param(
name="L_f",
init=L_f_init,
lo=L_f_bounds[0],
hi=L_f_bounds[1],
scale=L_f_init,
)
)
u_off_max = np.deg2rad(u_off_max_deg)
params.append(
Param(
name="u_off",
init=0.0,
lo=-u_off_max,
hi=u_off_max,
scale=float(np.deg2rad(1.0)),
)
)
if pin:
unknown = set(pin) - {p.name for p in params}
if unknown:
raise ValueError(
f"Cannot pin unknown parameters {sorted(unknown)}; this spec has "
f"{[p.name for p in params]}"
)
params = [
(p.model_copy(update={"init": pin[p.name], "pinned": True}) if p.name in pin else p)
for p in params
]
return ParamLayout(params=params)
[docs]
def predict_yaw(
w: WindowArrays,
x: NDArray[np.float64],
spec: FitSpec,
layout: ParamLayout,
irl_car: IrlCar,
y_init: float | None = None,
sign: float = 1.0,
) -> NDArray[np.float64]:
"""Predict one sub-window's yaw rate under ``spec`` at parameter vector ``x``.
The pipeline is: offset and mirror the input (``u_eff = sign * (u - u_off)``),
optionally lag it for front-tire relaxation, optionally decay each axle's
cornering stiffness against its own slip angle, then run the LPV biquad.
Parameters
----------
w : WindowArrays
The sub-window to predict over.
x : NDArray[np.float64]
Free-parameter vector, aligned with ``layout.free``.
spec : FitSpec
Selects the Ca scheduler and relaxation flavor.
layout : ParamLayout
Maps ``x`` (plus any pinned values) onto named parameters.
irl_car : IrlCar
Source of the mass, wheelbase, and CG split.
y_init : float, optional
Forced initial output; see ``simulate_lpv_yaw``.
sign : float, optional
Mirrors the input, for the sign-balanced fit views.
Returns
-------
NDArray[np.float64]
Predicted yaw rate (rad/s), same length as ``w.t``.
Notes
-----
The two optional extensions sit on top of the linear bicycle core (see
``bicycle_coeffs``). ``ca_scheduler="abs_alpha_per_axle"`` decays each
axle's cornering stiffness against its own measured slip-angle magnitude::
Ca_f(t) = Ca_f0 / (1 + alpha_f |alpha_f(t)|^p)
which is what lets one fit span the near-limit sub-windows where the tires
have left the linear region. ``relaxation="front"`` lags the steering input
by ``tau(V) = L_f / V``, modelling front-tire relaxation length; see
``apply_relaxation_lag``.
"""
chassis = ChassisParams(
Ca_f=layout.get(x, "Ca_f"),
Ca_r=layout.get(x, "Ca_r"),
Izz=layout.get(x, "Izz"),
)
u_eff = sign * (w.u_rad - layout.get(x, "u_off"))
if spec.relaxation == "front":
u_eff = apply_relaxation_lag(u_eff, w.speed, w.t, layout.get(x, "L_f"))
Ca_f_t = Ca_r_t = None
if spec.ca_scheduler != "none":
alpha_f_sig, alpha_r_sig = axle_slip_angles(w, irl_car)
p = layout.get(x, "p")
Ca_f_t = chassis.Ca_f / (1.0 + layout.get(x, "alpha_f") * np.power(alpha_f_sig, p))
Ca_r_t = chassis.Ca_r / (1.0 + layout.get(x, "alpha_r") * np.power(alpha_r_sig, p))
return simulate_lpv_yaw(
w.t, u_eff, w.speed, chassis, irl_car, Ca_f=Ca_f_t, Ca_r=Ca_r_t, y_init=y_init
)
[docs]
def build_fit_views(window_arrays: list[WindowArrays], mode: str) -> list[FitView]:
"""Build the ``(sub-window, sign, source id)`` views the residual is summed over.
Lapping data is rarely symmetric -- a track run one direction steers mostly
one way, which lets the optimizer trade a real steering offset against a
fake one. Mirroring views balances the excitation around zero.
Parameters
----------
window_arrays : list of WindowArrays
Sub-windows to build views from.
mode : {"none", "mirror_all", "flip_half"}
``"none"`` uses each sub-window once, unmirrored. ``"mirror_all"`` uses
each one twice, at both signs. ``"flip_half"`` mirrors the half of the
sub-windows with the most positive mean steer.
Returns
-------
list of FitView
Raises
------
ValueError
If ``mode`` is unrecognised.
"""
mode = mode.lower().strip()
if mode not in ("none", "mirror_all", "flip_half"):
raise ValueError(f"Unknown balance mode '{mode}'; expected none, mirror_all, or flip_half")
if mode == "none":
return [(w, +1.0, i) for i, w in enumerate(window_arrays)]
if mode == "mirror_all":
views: list[FitView] = []
for i, w in enumerate(window_arrays):
views.append((w, +1.0, i))
views.append((w, -1.0, i))
return views
mean_u = np.array([float(w.u_rad.mean()) for w in window_arrays])
most_positive = np.argsort(mean_u)[::-1][: len(window_arrays) // 2]
flipped = set(most_positive.tolist())
return [(w, -1.0 if i in flipped else +1.0, i) for i, w in enumerate(window_arrays)]
[docs]
def weighted_mean_u_deg(views: list[FitView]) -> float:
"""Sample-weighted mean steering angle across ``views`` (deg).
Near zero means the views are sign-balanced; far from zero means the
steering offset and the chassis parameters are fighting each other.
"""
num = 0.0
den = 0
for w, sign, _ in views:
num += float(sign * w.u_rad.mean() * len(w.u_rad))
den += len(w.u_rad)
if den == 0:
return float("nan")
return float(np.rad2deg(num / den))
def _window_metrics(w: WindowArrays, y_hat: NDArray[np.float64]) -> WindowMetrics:
"""Residual diagnostics for one sub-window's prediction."""
y_meas = w.yaw_rate
err_raw = y_meas - y_hat
dc = float(err_raw.mean())
err_ac = err_raw - dc
var_meas = float(np.var(y_meas))
return WindowMetrics(
subwindow_id=w.sw_id,
duration_s=w.duration_s,
mean_speed_mps=float(w.speed.mean()),
dc_offset_rad_s=dc,
rmse_raw_rad_s=float(np.sqrt(np.mean(err_raw**2))),
rmse_ac_rad_s=float(np.sqrt(np.mean(err_ac**2))),
yaw_rate_std_meas_rad_s=float(np.std(y_meas)),
vaf_ac_pct=float(100.0 * (1.0 - np.var(err_ac) / max(var_meas, 1e-12))),
)
def _predict_all(
window_arrays: list[WindowArrays],
x: NDArray[np.float64],
spec: FitSpec,
layout: ParamLayout,
irl_car: IrlCar,
n_params: int,
) -> Prediction:
"""Predict every sub-window at ``x`` and bundle the residual diagnostics.
Each sub-window is simulated in its own un-mirrored physical frame, so the
stored predictions always line up with the measured data regardless of how
the fit views were signed.
"""
y_hats: list[NDArray[np.float64]] = []
errs: list[NDArray[np.float64]] = []
rows: list[WindowMetrics] = []
n_samples = 0
rss = 0.0
for w in window_arrays:
y_meas = w.yaw_rate
y_hat = predict_yaw(w, x, spec, layout, irl_car, y_init=float(y_meas[0]))
err = y_meas - y_hat
y_hats.append(y_hat)
errs.append(err)
rows.append(_window_metrics(w, y_hat))
n_samples += len(y_meas)
rss += float(np.sum(err**2))
# AIC/BIC on a Gaussian likelihood, up to a constant common to every model.
mean_sq = max(rss, 1e-12) / max(n_samples, 1)
aic = 2.0 * n_params + n_samples * np.log(mean_sq)
bic = n_params * np.log(max(n_samples, 1)) + n_samples * np.log(mean_sq)
return Prediction(
y_hat_per_window=y_hats,
err_per_window=errs,
window_metrics=rows,
n_samples=n_samples,
aic=float(aic),
bic=float(bic),
)
[docs]
def fit_yaw(
fit_views: list[FitView],
window_arrays: list[WindowArrays],
spec: FitSpec,
layout: ParamLayout,
ctx: FitContext,
balance_mode: str = "none",
verbose: Literal[0, 1, 2] = 0,
max_nfev: int = 350,
) -> YawFit:
"""Fit ``spec`` to the sub-windows by bounded output-error least squares.
The optimizer minimises the concatenated residual across ``fit_views``.
Each view's simulation is initialised at that view's measured yaw rate, so
the residual starts at zero and the optimizer never has to fight an
initial-condition mismatch alongside the dynamics.
Parameters
----------
fit_views : list of FitView
Views the residual is summed over; see ``build_fit_views``.
window_arrays : list of WindowArrays
Source sub-windows, re-simulated once at the fitted parameters to
produce the reported predictions.
spec : FitSpec
The variant to fit.
layout : ParamLayout
The parameter vector and its bounds.
ctx : FitContext
Channel/geometry metadata, recorded on the fit.
balance_mode : str, optional
Not used here; only recorded on the fit so you can see later how the
views were weighted. See ``build_fit_views``.
verbose, max_nfev : optional
Forwarded to ``scipy.optimize.least_squares``.
Returns
-------
YawFit
The fitted parameters and their per-sub-window diagnostics.
"""
def residual(x: NDArray[np.float64]) -> NDArray[np.float64]:
chunks = []
for w, sign, _ in fit_views:
y_meas = sign * w.yaw_rate
y_hat = predict_yaw(w, x, spec, layout, ctx.irl_car, y_init=float(y_meas[0]), sign=sign)
chunks.append(y_meas - y_hat)
return np.concatenate(chunks)
result = least_squares(
residual,
layout.x0,
bounds=layout.bounds,
x_scale=layout.x_scale,
method="trf",
verbose=verbose,
max_nfev=max_nfev,
)
train = _predict_all(window_arrays, result.x, spec, layout, ctx.irl_car, layout.n_free)
return YawFit(
spec=spec,
layout=layout,
ctx=ctx,
x_hat=result.x,
train=train,
final_cost=float(result.cost),
n_eval=int(result.nfev),
status=int(result.status),
message=str(result.message),
balance_mode=balance_mode,
fit_view_count=len(fit_views),
source_window_count=len(window_arrays),
)
[docs]
def validate_fit(val_window_arrays: list[WindowArrays], fit: YawFit) -> Prediction:
"""Evaluate a frozen fit on held-out sub-windows -- no re-optimisation.
AIC / BIC use the fit's own parameter count, so they are directly
comparable against the training numbers on ``fit.train``.
Parameters
----------
val_window_arrays : list of WindowArrays
Held-out sub-windows, typically from another session's artifact.
fit : YawFit
The frozen fit to evaluate.
Returns
-------
Prediction
Predictions and metrics on ``val_window_arrays``.
"""
return _predict_all(
val_window_arrays,
fit.x_hat,
fit.spec,
fit.layout,
fit.ctx.irl_car,
fit.n_params,
)
[docs]
def compare_fits(fits: dict[str, YawFit]) -> FitComparison:
"""Table comparing several fits, keyed by display name."""
return FitComparison.from_fits(fits)
[docs]
def summarize_fit(fit: YawFit, sample_speeds: tuple[float, ...] = (10.0, 15.0, 20.0)) -> str:
"""Console summary of a fit plus its linear-region coefficients at ``sample_speeds``.
The coefficients are the *linear-region* ones, i.e. evaluated with the
Ca decay switched off, so they describe the car's small-slip behaviour.
"""
points = bicycle_points(fit.chassis_params(), fit.ctx.irl_car, sample_speeds)
lines = [str(fit), " Linear-region bicycle coefficients (Ca decay = 0):"]
lines.extend(f" {p}" for p in points)
return "\n".join(lines)