"""Understeer gradient measurement from steady-state test data."""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from suboptimumg.constants import G
from .derived import _require_col
if TYPE_CHECKING:
from .irl_car import IrlCar
[docs]
def measure_understeer_gradient(
df: pd.DataFrame,
time_range: tuple[float, float],
irl_car: IrlCar,
test_type: str = "constant_steer",
min_ay_g: float = 0.1,
force_zero_intercept: bool = False,
fit_mode: str = "linear",
symmetric: bool = False,
show_ramp: bool = False,
overlay_diagnostics: bool = False,
sw_col: str = "ludwig.steeringWheel.angle",
ay_col: str = "body.accLat",
kinematic_col: str = "body.steerAngle",
bicycle_col: str = "front.steerAngle.bicycle",
vel_col: str = "groundSpeed",
) -> go.Figure:
"""Measure the understeer gradient K from a steady-state test window.
Computes ``understeer_angle = |delta_bicycle| - |delta_kinematic|``
and fits a line against ``|a_y|`` in g's. The slope is the
understeer gradient K (deg/g). K > 0 = understeer, K < 0 = oversteer.
Works for both constant-speed/increasing-steer and
constant-steer/increasing-speed test procedures — the core math
is the same; *test_type* only affects the diagnostic subplot.
Parameters
----------
df : DataFrame
Must contain the columns referenced by the ``*_col`` arguments.
The DataFrame is **not** modified.
time_range : (t_start, t_end)
``time_s`` window to analyse.
irl_car : IrlCar
Used only for display / future extensions.
test_type : str
``"constant_steer"`` or ``"constant_speed"``. Controls which
variable is shown in the diagnostic subplot.
min_ay_g : float
Minimum lateral acceleration (g) to include in the fit.
Filters out low-speed noise.
force_zero_intercept : bool
If ``True``, force the fit through the origin (intercept = 0).
Useful when a systematic offset (e.g. steering sensor zero
error) is inflating the intercept.
fit_mode : str
``"linear"`` (default) — classic K = slope of UG vs a_y.
``"quadratic"`` — fits ``UG = c2 * a_y^2 + c1 * a_y + c0``.
``"cubic"`` — fits ``UG = c3 * a_y^3 + c2 * a_y^2 + c1 * a_y + c0``.
For non-linear modes the reported K is the local slope at
``a_y = 0`` (the ``c1`` term); higher-order terms capture
progressive under/oversteer with rising lateral load.
symmetric : bool
If ``True``, force the fitted equation to be symmetric across
the y-axis (``f(-x) = f(x)``) by zeroing all odd-degree terms.
Quadratic becomes ``c2 * a_y^2 + c0``; cubic likewise reduces
to even-only terms. Combined with ``force_zero_intercept``
the quadratic/cubic fits collapse to ``c2 * a_y^2``. Linear
+ symmetric is degenerate (constant only) and not allowed.
show_ramp : bool
If ``True``, add a subplot for the ramped variable (speed for
``constant_steer``, steering angle for ``constant_speed``).
overlay_diagnostics : bool
If ``True`` *and* ``show_ramp`` is ``True``, overlay both
diagnostic traces on a single subplot with dual y-axes instead
of separate subplots.
sw_col, ay_col, kinematic_col, bicycle_col, vel_col : str
Column names in *df*.
Returns
-------
go.Figure
Multi-panel figure: scatter + fit (top), diagnostic time-series
(bottom panel(s)).
"""
print(
"""
DEPRECATION WARNING: measure_understeer_gradient (loganalysis) is deprecated and will be removed in a future release.
Use suboptimumg.log_analysis.measure_understeer_gradient instead, which operates directly on a PERDA SingleRunData object.
"""
)
_require_col(df, bicycle_col, "Run compute_bicycle_steer_angle() first.")
_require_col(df, kinematic_col, "Run compute_kinematic_steer_angle() first.")
_require_col(df, ay_col, "Run compute_accelerations() first.")
_require_col(df, sw_col, "Steering wheel angle column not found.")
_require_col(df, vel_col, "Run compute_groundspeed() first.")
# --- Slice to window (local copy) ---
t_start, t_end = time_range
w = df.loc[t_start:t_end].copy()
if len(w) == 0:
raise ValueError(
f"No data in time_range=({t_start}, {t_end}). "
f"DataFrame index spans [{df.index.min():.1f}, {df.index.max():.1f}]."
)
# --- Compute understeer angle (absolute values for sign agnosticism) ---
delta_bicycle = np.abs(w[bicycle_col].values)
delta_kinematic = np.abs(w[kinematic_col].values)
understeer_deg = delta_bicycle - delta_kinematic
ay_g = np.abs(w[ay_col].values) / G
timestamps = w.index.values
# --- Filter out low-accel noise ---
mask = ay_g >= min_ay_g
ay_fit = ay_g[mask]
us_fit = understeer_deg[mask]
ts_fit = timestamps[mask]
if len(ay_fit) < 10:
raise ValueError(
f"Only {len(ay_fit)} points above min_ay_g={min_ay_g}g "
f"in window ({t_start}, {t_end}). Lower min_ay_g or widen the window."
)
# --- Fit (linear, quadratic, or cubic) ---
fit_mode = fit_mode.lower()
_fit_degrees = {"linear": 1, "quadratic": 2, "cubic": 3}
if fit_mode not in _fit_degrees:
raise ValueError(
f"fit_mode must be one of {list(_fit_degrees)}, got {fit_mode!r}"
)
deg = _fit_degrees[fit_mode]
if symmetric and fit_mode == "linear":
raise ValueError(
"symmetric=True with fit_mode='linear' is degenerate "
"(would reduce to a constant). Use quadratic or cubic."
)
# Build the set of polynomial powers to include in the fit.
# All powers from 0..deg, then drop odd powers if symmetric, then
# drop the constant term if force_zero_intercept.
powers = list(range(deg + 1))
if symmetric:
powers = [p for p in powers if p % 2 == 0]
if force_zero_intercept:
powers = [p for p in powers if p != 0]
if not powers:
raise ValueError(
"Fit constraints leave no free parameters "
"(force_zero_intercept + symmetric + linear?)."
)
n = len(ay_fit)
A = np.column_stack([ay_fit**p for p in powers])
sol, *_ = np.linalg.lstsq(A, us_fit, rcond=None)
# Reconstruct full coeffs vector for np.polyval (highest power first).
full = {p: 0.0 for p in range(deg + 1)}
for p, c in zip(powers, sol):
full[p] = float(c)
coeffs = np.array([full[p] for p in range(deg, -1, -1)])
c0 = full[0]
c1 = full.get(1, 0.0)
c2 = full.get(2, 0.0)
c3 = full.get(3, 0.0)
intercept = c0
K = c1 # local slope at a_y = 0
us_pred = np.polyval(coeffs, ay_fit)
ss_res = np.sum((us_fit - us_pred) ** 2)
ss_tot = np.sum((us_fit - np.mean(us_fit)) ** 2)
r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0
n_params = len(powers)
dof = n - n_params
std_err = (
np.sqrt(ss_res / dof) / np.sqrt(np.sum((ay_fit - ay_fit.mean()) ** 2))
if dof > 0
else float("inf")
)
intercept_str = " (forced)" if force_zero_intercept else f" = {intercept:.3f} deg"
extras = []
if 2 in full and (fit_mode != "linear"):
extras.append(f"c2 = {c2:.3f} deg/g²")
if 3 in full and fit_mode == "cubic":
extras.append(f"c3 = {c3:.3f} deg/g³")
extra_str = (", " + ", ".join(extras)) if extras else ""
sym_tag = " (symmetric)" if symmetric else ""
print(
f"[understeer] [{fit_mode}{sym_tag}] K = {K:.3f} deg/g{extra_str}, "
f"R² = {r_squared:.4f}, std_err = {std_err:.4f}, "
f"intercept{intercept_str}, n = {n}, "
f"window = {t_start:.1f}–{t_end:.1f} s"
)
# --- Equation title ---
_superscripts = {1: "", 2: "²", 3: "³"}
parts: list[str] = []
for p in sorted(powers, reverse=True):
c = full[p]
if p == 0:
term = f"{abs(c):.2f}"
elif p == 1:
term = f"{abs(c):.2f}x"
else:
term = f"{abs(c):.2f}x{_superscripts[p]}"
if not parts:
parts.append(("-" if c < 0 else "") + term)
else:
parts.append(("− " if c < 0 else "+ ") + term)
eq_title = "y = " + " ".join(parts)
time = w.index.values
sw_vals = w[sw_col].values
vel_vals = w[vel_col].values
if test_type == "constant_steer":
const_vals, const_label = sw_vals, "SW Angle (deg)"
ramp_vals, ramp_label = vel_vals, "Speed (m/s)"
else:
const_vals, const_label = vel_vals, "Speed (m/s)"
ramp_vals, ramp_label = sw_vals, "SW Angle (deg)"
# Determine layout
use_overlay = show_ramp and overlay_diagnostics
if use_overlay:
n_rows = 2
titles = [
f"Understeer Gradient — {eq_title} (R² = {r_squared:.3f})",
f"Diagnostics — {t_start:.0f}–{t_end:.0f} s",
]
specs = [[{}], [{"secondary_y": True}]]
heights = [0.60, 0.40]
elif show_ramp:
n_rows = 3
titles = [
f"Understeer Gradient — {eq_title} (R² = {r_squared:.3f})",
f"{const_label} (should be ~constant)",
f"{ramp_label} (ramped)",
]
specs = [[{}], [{}], [{}]]
heights = [0.50, 0.25, 0.25]
else:
n_rows = 2
titles = [
f"Understeer Gradient — {eq_title} (R² = {r_squared:.3f})",
f"{const_label} (should be ~constant) — {t_start:.0f}–{t_end:.0f} s",
]
specs = [[{}], [{}]]
heights = [0.65, 0.35]
fig = make_subplots(
rows=n_rows,
cols=1,
row_heights=heights,
vertical_spacing=0.12,
subplot_titles=titles,
specs=specs,
)
# Top panel: scatter + fit line
hover = [f"t = {t:.2f} s" for t in ts_fit]
fig.add_trace(
go.Scattergl(
x=ay_fit,
y=us_fit,
mode="markers",
marker=dict(size=3, opacity=0.4),
name="Data",
text=hover,
hovertemplate="a_y = %{x:.3f} g<br>UG = %{y:.2f} deg<br>%{text}<extra></extra>",
),
row=1,
col=1,
)
ay_line = np.linspace(ay_fit.min(), ay_fit.max(), 100)
fig.add_trace(
go.Scatter(
x=ay_line,
y=np.polyval(coeffs, ay_line),
mode="lines",
line=dict(color="red", width=2),
name=eq_title,
),
row=1,
col=1,
)
fig.update_xaxes(title_text="Lateral Acceleration (g)", row=1, col=1)
fig.update_yaxes(title_text="Understeer Angle (deg)", row=1, col=1)
# Diagnostic panel(s)
if use_overlay:
fig.add_trace(
go.Scattergl(
x=time, y=const_vals, mode="lines", name=const_label, line=dict(width=1)
),
row=2,
col=1,
secondary_y=False,
)
fig.add_trace(
go.Scattergl(
x=time, y=ramp_vals, mode="lines", name=ramp_label, line=dict(width=1)
),
row=2,
col=1,
secondary_y=True,
)
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_yaxes(title_text=const_label, row=2, col=1, secondary_y=False)
fig.update_yaxes(title_text=ramp_label, row=2, col=1, secondary_y=True)
elif show_ramp:
fig.add_trace(
go.Scattergl(
x=time, y=const_vals, mode="lines", name=const_label, line=dict(width=1)
),
row=2,
col=1,
)
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_yaxes(title_text=const_label, row=2, col=1)
fig.add_trace(
go.Scattergl(
x=time, y=ramp_vals, mode="lines", name=ramp_label, line=dict(width=1)
),
row=3,
col=1,
)
fig.update_xaxes(title_text="Time (s)", row=3, col=1)
fig.update_yaxes(title_text=ramp_label, row=3, col=1)
else:
fig.add_trace(
go.Scattergl(
x=time, y=const_vals, mode="lines", name=const_label, line=dict(width=1)
),
row=2,
col=1,
)
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_yaxes(title_text=const_label, row=2, col=1)
total_height = 700 if n_rows <= 2 else 850
fig.update_layout(
height=total_height,
showlegend=True,
template="plotly_white",
)
return fig