[education/kstars] kstars/ekos/guide/offlinetrainer: offline_trainer: extract shared pulse-response fitting, add PID auto-tune and QC tooling
Jasem Mutlaq <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit 3ba211d9837fc74e75e7cc724911143cf58699df by Jasem Mutlaq. Committed on 04/08/2026 at 02:59. Pushed by mutlaqja into branch 'master'. offline_trainer: extract shared pulse-response fitting, add PID auto-tune and QC tooling - New pulse_response_fit.py: extracts the step-response curve-fitting logic (previously private to train_harmonic.py's kappa/tau spring fit) into a shared, mount-agnostic fit_pulse_response(), now also consumed by pid_autotune.py's SIMC-style gain recommendation. - New pid_autotune.py: recommend_pid_gains(), a FOPDT/SIMC-style base RA/DEC proportional(+integral) gain recommendation from pulse-response sessions, for any mount type -- the Python source of truth behind the live C++ gain-lock (AIGuideProtocol::applyPIDAutoTuneGainLock()). - New validate_sysid.py: fast QC pass over a sysid_data_*.json file (works mid-run or finished), catching data-quality problems (all-zero frames, insufficient fits, coverage gaps) before a full training run. - New analyze_oscillation.py / tail_ai_log.py: post-hoc oscillation diagnosis and a live tailer/plotter for the AI debug CSV. - train_harmonic.py: the local kappa/tau spring-fit implementation moves to pulse_response_fit.py; the fit call itself is commented out and replaced with hardcoded defaults (kappa=0.0, tau=1.5) since it has never resolved above the noise floor on any rig tested -- see the physical analysis in offlinetrainer/AI_ARCHITECTURE.md §6.2.1. Also wires in pid_autotune.recommend_pid_gains() for the base-gain cross-check. - train.py / train_direct_drive.py / train_worm_gear.py / README.md: minor updates to match. Co-Authored-By: Claude Sonnet 5 <[email protected]> M +15 -7 kstars/ekos/guide/offlinetrainer/README.md A +393 -0 kstars/ekos/guide/offlinetrainer/analyze_oscillation.py A +190 -0 kstars/ekos/guide/offlinetrainer/pid_autotune.py A +243 -0 kstars/ekos/guide/offlinetrainer/pulse_response_fit.py A +252 -0 kstars/ekos/guide/offlinetrainer/tail_ai_log.py M +12 -9 kstars/ekos/guide/offlinetrainer/train.py M +27 -1 kstars/ekos/guide/offlinetrainer/train_direct_drive.py M +35 -387 kstars/ekos/guide/offlinetrainer/train_harmonic.py M +27 -1 kstars/ekos/guide/offlinetrainer/train_worm_gear.py A +336 -0 kstars/ekos/guide/offlinetrainer/validate_sysid.py https://invent.kde.org/education/kstars/-/commit/3ba211d9837fc74e75e7cc724911143cf58699df diff --git a/kstars/ekos/guide/offlinetrainer/README.md b/kstars/ekos/guide/offlinetrainer/README.md index edb9061d98..1ec6aedb45 100644 --- a/kstars/ekos/guide/offlinetrainer/README.md +++ b/kstars/ekos/guide/offlinetrainer/README.md @@ -32,7 +32,7 @@ python train.py --sysid-data ./sysid_data.json --output ./weights.json * `--mount-type <type>`: Force a specific mount type (`WORM_GEAR`, `HARMONIC_DRIVE`, `DIRECT_DRIVE`). Usually not needed as it auto-detects from the sysid file. * `--verbose`: Print detailed logs of the physics fitting and the MLP training epochs. * `--gpu`: Use CUDA/MPS to accelerate training if available (not necessary for these small models). -* `--pid-lambda-factor <float>`: HARMONIC_DRIVE only. Design parameter for the advisory PID auto-tune recommendation below (default `3.0`); larger values recommend a slower, more conservative gain. +* `--pid-lambda-factor <float>`: Design parameter for the advisory PID auto-tune recommendation below (default `3.0`); larger values recommend a slower, more conservative gain. Applies to any mount type whose sysid data includes `pulse_response` sessions. ## 4. Loading the Model into KStars @@ -41,17 +41,25 @@ python train.py --sysid-data ./sysid_data.json --output ./weights.json 3. Under the **AI Guiding** section, set the **Weights File** path to point to your `weights.json`. 4. Start guiding. The AI will automatically engage! -## 5. PID Auto-Tune Recommendation (HARMONIC_DRIVE, advisory only) +## 5. PID Auto-Tune Recommendation (all mount types, advisory only) -For `HARMONIC_DRIVE` mounts, the trainer also derives a recommended base RA/DEC -proportional (+ conservative integral) gain from the same `pulse_response` -step-response sessions used to fit κ/τ, via a conservative SIMC/IMC-style -step-response tuning rule. This is printed to the console and saved into -`weights.json` as `recommended_ra_proportional_gain`, `recommended_ra_integral_gain`, +For any mount type whose sysid data includes `pulse_response` sessions, the +trainer also derives a recommended base RA/DEC proportional (+ conservative +integral) gain from that step-response data, via a conservative SIMC/IMC-style +step-response tuning rule (`pid_autotune.py`, shared by all three trainers). +This is printed to the console and saved into `weights.json` as +`recommended_ra_proportional_gain`, `recommended_ra_integral_gain`, `recommended_dec_proportional_gain`, `recommended_dec_integral_gain`, plus a `pid_autotune` block with the underlying fit diagnostics (process gain, τ, dead time, confidence). +`pulse_response` sessions are collected by a protocol phase (**PID Auto-Tune +(Step-Response Gain Calibration)** in AI Guiding settings, enabled by default — +it adds roughly 10 minutes to the AI Guide Wizard run). If disabled, +`recommended_*` fields are `null` and `pid_autotune` reports +`"confidence": "unavailable"` per axis — this is the normal, expected result +for a run collected with the option off. + **This is advisory only and is never applied automatically.** Nothing in KStars reads these fields back — review the numbers (especially the `confidence` flag: `"unavailable"` means don't use them, `"low"` means the dead time could not be diff --git a/kstars/ekos/guide/offlinetrainer/analyze_oscillation.py b/kstars/ekos/guide/offlinetrainer/analyze_oscillation.py new file mode 100644 index 0000000000..2b0ce3ed0d --- /dev/null +++ b/kstars/ekos/guide/offlinetrainer/analyze_oscillation.py @@ -0,0 +1,393 @@ +#!/usr/bin/env python3 +""" +analyze_oscillation.py — Oscillation diagnostics for ACTIVE-mode AI Guider sessions. + +Unlike evaluate_shadow.py (which only handles SHADOW-mode counterfactual data), this +script reads the AI debug CSV from a real closed-loop session — AI-active, Shadow-mode +baseline, or plain standard guiding with Shadow Mode logging enabled — and helps answer: + + Is the AI+PID blend double-correcting (loop-gain / resonance), or is it a PE + phase mismatch (the AI confidently pushing at the wrong point in the cycle)? + +It relies on the blended-pulse breakdown columns (ra_algorithm, ra_prop_response_ms, +ra_integral_response_ms, ra_ai_response_ms, ra_active_prop_gain, ra_total_pulse_ms, +ra_direction, ra_suppressed, and the dec_* equivalents) that record what was ACTUALLY +sent to the mount each frame, not just the AI's internal prediction. Older CSVs (from +before this logging was added) lack these columns; the script degrades gracefully and +skips the blend-specific analyses for them. + +Frame quality policy is adapted from EVALUATION.md Sec 4a in the ekos-ai-guider design +repo (dropped frames via dt outliers, restart-gap re-convergence, error spikes, saturated +pulses). Two columns that policy also uses — ErrorCode and SNR — are not present in this +CSV (they live in the separate PHD2-format guide log), so those two criteria are omitted +here; note this if you need bullet-proof frame exclusion. + +The ringing diagnostic is adapted from acf_plots.py's autocorrelation approach (also in +that design repo) rather than an ad hoc sign-flip counter: a genuinely oscillating +closed loop shows strong NEGATIVE autocorrelation at lag 1 (and often lag 2) in the +actual pulse series — a cleaner statistical signature than counting sign flips. + +Usage: + python3 analyze_oscillation.py --log path/to/ai_guider_*.csv + python3 analyze_oscillation.py --logdir ~/.local/share/kstars/ai_debug_logs/ --latest 1 + python3 analyze_oscillation.py --log active.csv --baseline standard_baseline.csv + python3 analyze_oscillation.py --log active.csv --weights weights_harmonic.json +""" + +import argparse +import glob +import json +import os +import sys + +import numpy as np +import pandas as pd + +BLEND_COLUMNS = [ + "ra_algorithm", "ra_prop_response_ms", "ra_integral_response_ms", "ra_ai_response_ms", + "ra_active_prop_gain", "ra_total_pulse_ms", "ra_direction", "ra_suppressed", + "dec_algorithm", "dec_prop_response_ms", "dec_integral_response_ms", "dec_ai_response_ms", + "dec_active_prop_gain", "dec_total_pulse_ms", "dec_direction", "dec_suppressed", +] + +REQUIRED_COLUMNS = [ + "t_session", "dt", "ra_error_arcsec", "dec_error_arcsec", "conf", "ai_state", +] + + +def load_csv(path: str) -> pd.DataFrame: + """Load an AI debug CSV. Tags df.attrs['has_blend_columns'] for graceful degradation.""" + df = pd.read_csv(path) + + missing = [c for c in REQUIRED_COLUMNS if c not in df.columns] + if missing: + raise ValueError(f"CSV missing required columns: {missing}") + + has_blend = all(c in df.columns for c in BLEND_COLUMNS) + df.attrs["has_blend_columns"] = has_blend + if not has_blend: + print(f" [warn] {os.path.basename(path)}: old-format CSV (no blend breakdown columns) — " + "skipping blend-specific analyses (algorithm/total-pulse breakdown, ACF on pulse series).", + file=sys.stderr) + return df + + +# ── Frame quality policy (adapted from EVALUATION.md Sec 4a) ───────────────────────────── +def classify_frames(df: pd.DataFrame) -> pd.Series: + """EXCLUDED mask: True = exclude. Dropped frames (dt outlier) + post-restart re-convergence. + + Note: ErrorCode and SNR criteria from EVALUATION.md Sec 4a are not applicable here — this + CSV does not carry those columns (they live in the separate PHD2-format guide log). + """ + median_dt = df["dt"].median() + excluded = df["dt"] > 3 * median_dt + + # Restart gap: dt > 30s marks a guide abort/restart; the first 10 frames after are + # re-convergence, not steady-state behavior. + gap_idx = df.index[df["dt"] > 30.0] + restart_frames = set() + for idx in gap_idx: + pos = df.index.get_loc(idx) + for i in range(pos, min(pos + 10, len(df))): + restart_frames.add(df.index[i]) + excluded = excluded | df.index.isin(restart_frames) + + return excluded + + +def classify_flagged(df: pd.DataFrame, valid_mask: pd.Series) -> pd.Series: + """Among valid frames: error spikes (excluded from ringing/periodogram analysis only — + they're real events but break the linear/steady-state assumptions). + + Deliberately NOT flagging "saturated pulses" here the way EVALUATION.md Sec 4a does: + that policy uses pulse_ms >= 0.95 * the CONFIGURED max-pulse ceiling, which isn't logged + per-frame in this CSV. Using the session's own observed max pulse as a proxy instead + would backfire on exactly the failure mode this script hunts for — a session genuinely + oscillating at sustained high amplitude would have most of its frames sit near its own + max and get flagged out, hiding the ringing signal. So only the error-spike criterion + (scaled off the error series, not the pulse series) is applied here. + """ + valid = df[valid_mask] + if len(valid) == 0: + return pd.Series(False, index=df.index) + + ra_rms = float(np.sqrt(np.mean(valid["ra_error_arcsec"] ** 2))) + dec_rms = float(np.sqrt(np.mean(valid["dec_error_arcsec"] ** 2))) + + flagged = (df["ra_error_arcsec"].abs() > 5 * max(ra_rms, 1e-6)) | \ + (df["dec_error_arcsec"].abs() > 5 * max(dec_rms, 1e-6)) + + return valid_mask & flagged + + +# ── Confidence segmentation ─────────────────────────────────────────────────────────────── +def segment_by_confidence(df: pd.DataFrame, valid_mask: pd.Series, min_confidence: float = 0.5) -> dict: + valid = df[valid_mask] + if len(valid) == 0: + return {"error": "No valid frames after quality filtering"} + + confident = valid[(valid["ai_state"].isin(["ACTIVE", "SHADOW"])) & (valid["conf"] >= min_confidence)] + not_confident = valid[~valid.index.isin(confident.index)] + + def axis_rms(frame, col): + return float(np.sqrt(np.mean(frame[col] ** 2))) if len(frame) else float("nan") + + return { + "total_frames": len(df), + "valid_frames": len(valid), + "confident_frames": len(confident), + "not_confident_frames": len(not_confident), + "ra_rms_confident": axis_rms(confident, "ra_error_arcsec"), + "dec_rms_confident": axis_rms(confident, "dec_error_arcsec"), + "ra_rms_not_confident": axis_rms(not_confident, "ra_error_arcsec"), + "dec_rms_not_confident": axis_rms(not_confident, "dec_error_arcsec"), + "mean_confidence": float(confident["conf"].mean()) if len(confident) else float("nan"), + } + + +# ── Periodogram (Lomb-Scargle, same tool train_harmonic.py uses for PE detection) ──────── +def periodogram_report(df: pd.DataFrame, valid_mask: pd.Series, weights_path: str = None) -> dict: + valid = df[valid_mask].sort_values("t_session") + if len(valid) < 20: + return {"error": "Too few valid frames for a periodogram"} + + try: + import scipy.signal + import scipy.stats + except ImportError: + return {"error": "scipy is required for periodogram_report (pip install scipy)"} + + t = valid["t_session"].values.astype(float) + span = t[-1] - t[0] + if span <= 0: + return {"error": "Degenerate time span"} + + dt_med = float(np.median(np.diff(t))) if len(t) > 1 else 2.0 + nyquist = 0.5 / max(dt_med, 1e-3) + f_min = max(2.0 / span, 0.002) + f_max = min(0.5, nyquist * 0.9) + if f_min >= f_max: + return {"error": "Session too short to resolve any periodic structure"} + + f_search = np.geomspace(f_min, f_max, 4000) + omega = 2 * np.pi * f_search + + result = {} + for axis, col in (("ra", "ra_error_arcsec"), ("dec", "dec_error_arcsec")): + y = valid[col].values.astype(float) + slope, intercept, _, _, _ = scipy.stats.linregress(t, y) + y_detrended = y - (slope * t + intercept) + power = scipy.signal.lombscargle(t, y_detrended, omega, precenter=True) + noise_floor = np.median(power) + 1e-12 + + top_idx = np.argsort(power)[-3:][::-1] + top_peaks = [(float(1.0 / f_search[i]), float(power[i] / noise_floor)) for i in top_idx] + result[axis] = {"top_peaks_period_s_snr": top_peaks} + + if weights_path and os.path.exists(weights_path): + try: + with open(weights_path) as f: + weights = json.load(f) + if weights.get("mount_type") == "HARMONIC_DRIVE": + phys = weights.get("physical", {}) + pe_period = phys.get("pe_period", 0.0) + pe2_period = phys.get("pe2_period", 0.0) + result["trained_pe_periods_s"] = [p for p in (pe_period, pe2_period) if p and p > 0] + except (json.JSONDecodeError, OSError) as e: + print(f" [warn] could not read --weights {weights_path}: {e}", file=sys.stderr) + + return result + + +# ── ACF-based ringing diagnostic (adapted from acf_plots.py) ───────────────────────────── +def compute_autocorr(series: np.ndarray, max_lag: int) -> np.ndarray: + n = len(series) + s = series - series.mean() + ac = np.correlate(s, s, mode="full") + ac = ac[n - 1: n - 1 + max_lag] + denom = ac[0] if ac[0] != 0 else 1e-12 + return ac / denom + + +def acf_ringing_report(df: pd.DataFrame, valid_mask: pd.Series, max_lag: int = 10) -> dict: + """Strong NEGATIVE autocorrelation at lag 1 (sign flipping every frame) or lag 2 is the + statistical signature of closed-loop ringing/double-correction. Applied to the residual + (ra/dec_error_arcsec) and — when available — the actual total pulse sent, which is the + more direct "is the mount being oscillated" signal. + """ + valid = df[valid_mask] + if len(valid) < max_lag + 5: + return {"error": "Too few valid frames for autocorrelation"} + + n = len(valid) + sig = 1.96 / np.sqrt(n) + + report = {"n_frames": n, "significance_bound": float(sig), "series": {}} + + series_map = {"ra_error": valid["ra_error_arcsec"].values, "dec_error": valid["dec_error_arcsec"].values} + if df.attrs.get("has_blend_columns"): + series_map["ra_total_pulse"] = valid["ra_total_pulse_ms"].values + series_map["dec_total_pulse"] = valid["dec_total_pulse_ms"].values + + for name, series in series_map.items(): + ac = compute_autocorr(np.asarray(series, dtype=float), max_lag) + lag1, lag2 = float(ac[1]), float(ac[2]) if max_lag > 2 else float("nan") + ringing = bool(lag1 < -sig or (max_lag > 2 and lag2 < -sig)) + report["series"][name] = { + "lag1": lag1, + "lag2": lag2, + "ringing_signature": ringing, + } + + return report + + +# ── Session comparison ──────────────────────────────────────────────────────────────────── +def compare_sessions(primary: dict, baseline: dict) -> dict: + def improvement(base, cand): + if base is None or cand is None or base != base or cand != cand or base < 1e-9: + return float("nan") + return (1.0 - cand / base) * 100.0 + + return { + "ra_improvement_pct": improvement(baseline.get("ra_rms_confident"), primary.get("ra_rms_confident")), + "dec_improvement_pct": improvement(baseline.get("dec_rms_confident"), primary.get("dec_rms_confident")), + } + + +# ── Reporting ────────────────────────────────────────────────────────────────────────────── +def print_report(label: str, seg: dict, periodogram: dict, acf: dict): + print("=" * 78) + print(f" AI GUIDER — OSCILLATION DIAGNOSTIC: {label}") + print("=" * 78) + + if "error" in seg: + print(f" {seg['error']}") + return + print(f" Total frames : {seg['total_frames']}") + print(f" Valid frames : {seg['valid_frames']}") + print(f" Confident frames : {seg['confident_frames']} (mean conf {seg['mean_confidence']:.2f})") + print(f" Not-confident frames: {seg['not_confident_frames']}") + print() + print(f" {'Axis':<6} {'RMS (confident)':>18} {'RMS (not confident)':>22}") + print(f" {'RA':<6} {seg['ra_rms_confident']:>17.3f}\" {seg['ra_rms_not_confident']:>21.3f}\"") + print(f" {'DEC':<6} {seg['dec_rms_confident']:>17.3f}\" {seg['dec_rms_not_confident']:>21.3f}\"") + print() + + if "error" in periodogram: + print(f" Periodogram: {periodogram['error']}") + else: + trained = periodogram.get("trained_pe_periods_s") + for axis in ("ra", "dec"): + peaks = periodogram[axis]["top_peaks_period_s_snr"] + desc = ", ".join(f"{p:.1f}s (SNR {s:.0f})" for p, s in peaks) + print(f" {axis.upper()} residual periodogram top peaks: {desc}") + if trained: + print(f" Trained PE period(s): {', '.join(f'{p:.1f}s' for p in trained)}") + print(" -> Power concentrated near a trained PE period suggests PHASE MISMATCH.") + print(" -> Power concentrated elsewhere (short/near-Nyquist) suggests LOOP-GAIN /") + print(" double-correction (AI + P fighting each other), not a PE tracking issue.") + print() + + if "error" in acf: + print(f" ACF ringing check: {acf['error']}") + else: + print(f" ACF ringing check (significance bound ±{acf['significance_bound']:.3f}):") + for name, s in acf["series"].items(): + flag = " <-- RINGING SIGNATURE" if s["ringing_signature"] else "" + print(f" {name:<16} lag1={s['lag1']:+.3f} lag2={s['lag2']:+.3f}{flag}") + any_ringing = any(s["ringing_signature"] for s in acf["series"].values()) + pulse_ringing = any(s["ringing_signature"] for n, s in acf["series"].items() if "pulse" in n) + error_ringing = any(s["ringing_signature"] for n, s in acf["series"].items() if "error" in n) + print() + if pulse_ringing: + print(" VERDICT: the actual pulse sent to the mount is ringing (strong negative lag-1/2") + print(" autocorrelation) — consistent with AI+PID double-correction. Check") + print(" Options::aIProportionalBackoff / aIPredictionGain, or try the AI-only") + print(" reload with a lower aIPredictionGain via reloadAIWeights() and re-watch.") + elif error_ringing: + print(" VERDICT: the residual error is ringing but the pulse series isn't (or isn't") + print(" logged) — check the PE-period periodogram above for a phase-mismatch signature.") + elif any_ringing: + print(" VERDICT: some ringing signature detected — see per-series breakdown above.") + else: + print(" VERDICT: no autocorrelation ringing signature detected in this session.") + print() + print("-" * 78) + print() + + +def analyze_one(path: str, weights_path: str, min_confidence: float, max_lag: int): + df = load_csv(path) + excluded = classify_frames(df) + valid_mask = ~excluded + flagged = classify_flagged(df, valid_mask) + analysis_mask = valid_mask & ~flagged + + seg = segment_by_confidence(df, valid_mask, min_confidence) + periodogram = periodogram_report(df, analysis_mask, weights_path) + acf = acf_ringing_report(df, analysis_mask, max_lag) + return seg, periodogram, acf + + +def main(): + parser = argparse.ArgumentParser( + description="Diagnose AI Guider oscillations from ACTIVE-mode (or Shadow-mode baseline) debug CSVs.") + parser.add_argument("--log", type=str, nargs="*", help="Path(s) to AI debug CSV files (glob supported)") + parser.add_argument("--logdir", type=str, default=None, help="Directory containing AI debug CSVs") + parser.add_argument("--latest", type=int, default=0, help="Only analyze the N most recent CSVs (0=all)") + parser.add_argument("--baseline", type=str, default=None, + help="A second CSV (e.g. Shadow-mode/standard-guiding session) for side-by-side RMS") + parser.add_argument("--weights", type=str, default=None, + help="weights_*.json to annotate the periodogram with the trained PE period(s)") + parser.add_argument("--min-confidence", type=float, default=0.5, + help="Confidence threshold splitting AI-active vs not-yet-confident frames (default 0.5)") + parser.add_argument("--flip-window", type=int, default=5, + help="Max lag (frames) for the ACF ringing check (default 5)") + args = parser.parse_args() + + csv_files = [] + if args.log: + for pattern in args.log: + csv_files.extend(glob.glob(os.path.expanduser(pattern))) + if args.logdir: + logdir = os.path.expanduser(args.logdir) + csv_files.extend(glob.glob(os.path.join(logdir, "ai_guider_*.csv"))) + if not csv_files: + default_dir = os.path.expanduser("~/.local/share/kstars/ai_debug_logs/") + csv_files = glob.glob(os.path.join(default_dir, "ai_guider_*.csv")) + if not csv_files: + print("No AI debug CSV files found. Specify --log or --logdir.", file=sys.stderr) + sys.exit(1) + + csv_files = sorted(set(csv_files), key=os.path.getmtime) + if args.latest > 0: + csv_files = csv_files[-args.latest:] + + print(f"Found {len(csv_files)} CSV file(s) to analyze.") + print() + + baseline_seg = None + if args.baseline: + try: + baseline_seg, _, _ = analyze_one(args.baseline, args.weights, args.min_confidence, args.flip_window) + print_report(f"BASELINE: {os.path.basename(args.baseline)}", baseline_seg, {"error": "n/a for baseline"}, + {"error": "n/a for baseline"}) + except Exception as e: + print(f"Could not analyze baseline {args.baseline}: {e}", file=sys.stderr) + + for path in csv_files: + try: + seg, periodogram, acf = analyze_one(path, args.weights, args.min_confidence, args.flip_window) + print_report(os.path.basename(path), seg, periodogram, acf) + if baseline_seg and "error" not in seg: + cmp = compare_sessions(seg, baseline_seg) + print(f" vs. baseline: RA {cmp['ra_improvement_pct']:+.1f}% " + f"DEC {cmp['dec_improvement_pct']:+.1f}%") + print() + except Exception as e: + print(f"[{os.path.basename(path)}] ERROR: {e}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/kstars/ekos/guide/offlinetrainer/pid_autotune.py b/kstars/ekos/guide/offlinetrainer/pid_autotune.py new file mode 100644 index 0000000000..12c66d7e79 --- /dev/null +++ b/kstars/ekos/guide/offlinetrainer/pid_autotune.py @@ -0,0 +1,190 @@ +""" +offline_trainer/pid_autotune.py — recommend base RA/DEC PID guide gains from +pulse_response step-response data, for any mount type. + +Offline, calibration-time-only calculation (pid_autotune_plan.md). Reuses the +per-fit P_fit/tau_fit/residual_std/t_arr values pulse_response_fit.fit_pulse_ +response() already computes from pulse_response sessions -- no new data +collection, just a second consumer of the same step-response fits. Extracts a +per-axis FOPDT-like model (process gain K, time constant tau, dead time L), +applies a conservative SIMC/IMC-style PI tuning rule, and back-converts the +result into KStars' dimensionless proportional_gain convention via the +calibrated ms_per_arcsec rate recorded in standard_guiding sessions. + +This is advisory only. Nothing in KStars or any trainer applies the result +automatically -- callers surface it as recommended_ra/dec_proportional_gain +(+ integral_gain) fields plus a confidence flag, for a human to review before +manually updating Options::rA/dECProportionalGain(). + +Originally implemented inside train_harmonic.py (the only mount type with +pulse_response data at the time); moved here, unchanged, once WORM_GEAR and +DIRECT_DRIVE gained their own pulse-response phases -- the calculation itself +was never Harmonic-Drive-specific, only the data collection was. + +SPDX-License-Identifier: GPL-2.0-or-later +""" + +import numpy as np + +from pulse_response_fit import fit_pulse_response + +# Default SIMC design parameter: lambda = max(tau, SIMC_LAMBDA_L_FACTOR * dead_time). +# Larger -> slower/more robust closed loop; smaller -> faster/less margin. Kept as +# an explicit, overridable constant per the plan's "expose lambda as a documented +# knob, not a hidden constant" (pid_autotune_plan.md §3.2). +SIMC_LAMBDA_L_FACTOR = 3.0 + +# KStars' "integral gain" (Options::rAIntegralGain()/dECIntegralGain()) multiplies +# a ~100-frame moving average of drift (gmath.cpp::processAxis(), drift_integral[k]), +# not a classical accumulating integrator. SIMC's tau_I (a reset *time*) has no +# principled mapping onto that EMA-style term, so rather than inventing one, the +# integral recommendation is a conservative fixed fraction of the proportional +# recommendation (pid_autotune_plan.md §1: "keep I conservative"). tau_I is still +# reported in the output for reference. +INTEGRAL_GAIN_CONSERVATIVE_FRACTION = 0.25 + +# Below this many usable step-response fits, the recommendation is flagged "low" +# confidence rather than withheld outright -- still informative, not to be applied +# unattended (pid_autotune_plan.md §4). +MIN_FITS_FOR_MEDIUM_CONFIDENCE = 6 + + +def _effective_pixel_scale(sysid): + """Pixel scale in arcsec/px. Older exports recorded it without the binning factor.""" + eq = sysid.get("equipment", {}) + ps = float(eq.get("pixel_scale_arcsec_per_px", 1.0) or 1.0) + if not eq.get("pixel_scale_includes_binning", False): + b = str(sysid.get("model_fingerprint", {}).get("guide_binning", "1x1")) + try: + bf = max(1, int(b.split("x")[0])) + except ValueError: + bf = 1 + if bf > 1: + ps *= bf + return ps + + +def _estimate_dead_time_s(t_arr: np.ndarray, pos_arr: np.ndarray, residual_std: float) -> float: + """ + First t at which |pos(t)| exceeds ~2.5x the fit's residual noise -- a proxy + for FOPDT dead time L (pid_autotune_plan.md §3.1). If even the first sample + already exceeds threshold (the common case at ~2-3s/frame cadence), this + returns t_arr[0]: L is only known to be <= the first sample, not resolved any + finer -- see _recommend_axis_pid_gain()'s "resolution_limited" flag. + """ + threshold = max(2.5 * residual_std, 1e-6) + for t, p in zip(t_arr, pos_arr): + if abs(p) > threshold: + return float(t) + return float(t_arr[-1]) + + +def _calibration_ms_per_arcsec(sysid: dict, axis: str) -> float: + """ + Median calibrated ms-per-arcsec pulse rate for this axis, from whichever + sessions recorded it (standard_guiding sessions carry ra_ms_per_arcsec/ + dec_ms_per_arcsec). This is Calibration::ra/decPulseMillisecondsPerArcsecond() + at collection time -- the same normalization gmath.cpp::processAxis() applies + to proportional_gain, needed to convert a physical SIMC Kc (ms/arcsec) back + into KStars' dimensionless aggressiveness (pid_autotune_plan.md §3.3). + """ + key = "ra_ms_per_arcsec" if axis.upper() == "RA" else "dec_ms_per_arcsec" + values = [s[key] for s in sysid["sessions"] if s.get(key, 0.0) and s[key] > 0.0] + return float(np.median(values)) if values else 0.0 + + +def _recommend_axis_pid_gain(sysid: dict, axis: str, guide_exp: float, + lambda_l_factor: float, verbose: bool) -> dict: + """ + Derive a recommended base proportional/integral gain for one axis from its + pulse_response step-response fits, via a conservative SIMC/IMC-style PI rule + (pid_autotune_plan.md §3.2). Mount-agnostic: works for any mount type whose + sysid data includes pulse_response sessions for this axis. Returns a dict; + confidence "unavailable" means the numbers (if present at all) should not + be used -- most commonly because this mount type/run has no pulse_response + data yet (Options::aIPIDAutoTune() was off, or this mount class's + protocol doesn't collect it). + """ + _, _, fit_info = fit_pulse_response(sysid, axis, guide_exp, False, return_fits=True) + fits = fit_info["fits"] + + if fit_info["sign_consistent"] is False: + return {"confidence": "unavailable", + "reason": "pulse-response signs inconsistent across directions -- fits are noise, not mechanics"} + if not fits: + return {"confidence": "unavailable", "reason": "no usable pulse-response step-response fits"} + + cal_ms_per_arcsec = _calibration_ms_per_arcsec(sysid, axis) + if cal_ms_per_arcsec <= 0.0: + return {"confidence": "unavailable", + "reason": "no calibrated ms_per_arcsec recorded for this axis (need a standard_guiding session)"} + + pixel_scale = _effective_pixel_scale(sysid) + + K_samples, tau_samples, L_samples, t_first_samples = [], [], [], [] + for f in fits: + if f["pulse_magnitude_ms"] <= 0.0: + continue + K_samples.append(abs(f["P_fit_px"]) * pixel_scale / f["pulse_magnitude_ms"]) + t_first_samples.append(f["t_first_s"]) + L_samples.append(_estimate_dead_time_s(f["t_arr"], f["pos_arr"], f["residual_std_px"])) + if f["tau_fit_s"] <= 9.8: # same "pinned at bound, degenerate with drift" guard as the spring fit + tau_samples.append(f["tau_fit_s"]) + + if not K_samples: + return {"confidence": "unavailable", "reason": "no fit had a usable pulse magnitude"} + + K = float(np.median(K_samples)) # arcsec of steady-state response per ms of pulse + L = float(np.median(L_samples)) + # tau can't be resolved below the sampling floor either: a fit whose tau_fit + # landed below t_first ("spring already released") is only known to be + # <= t_first, not physically ~0 -- floor it the same way L is floored, so a + # spuriously tiny fitted tau can't produce a dangerously aggressive Kc. + tau_raw = float(np.median(tau_samples)) if tau_samples else float(np.median(t_first_samples)) + tau = max(tau_raw, L) + resolution_limited = bool(np.isclose(L, float(np.median(t_first_samples)), rtol=0.05)) + + lam = max(tau, lambda_l_factor * L) + Kc = (1.0 / K) * tau / (lam + L) # ms of pulse per arcsec of error + tau_i = min(tau, 4.0 * (lam + L)) # SIMC reset time, reported only -- see constant doc above + + proportional_gain = Kc / cal_ms_per_arcsec + integral_gain = INTEGRAL_GAIN_CONSERVATIVE_FRACTION * proportional_gain + confidence = "low" if (resolution_limited or len(fits) < MIN_FITS_FOR_MEDIUM_CONFIDENCE) else "medium" + + if verbose: + print(f" [{axis} PID] K={K:.5f} arcsec/ms tau={tau:.2f}s L={L:.2f}s lambda={lam:.2f}s " + f"(n={len(fits)} fits, cal={cal_ms_per_arcsec:.1f}ms/arcsec)") + print(f" [{axis} PID] Recommended proportional_gain={proportional_gain:.3f} " + f"integral_gain={integral_gain:.3f} confidence={confidence}") + + return { + "confidence": confidence, + "proportional_gain": float(np.clip(proportional_gain, 0.0, 1.0)), + "integral_gain": float(np.clip(integral_gain, 0.0, 1.0)), + "process_gain_arcsec_per_ms": K, + "tau_s": tau, + "dead_time_s": L, + "lambda_s": lam, + "tau_i_s": tau_i, + "calibration_ms_per_arcsec": cal_ms_per_arcsec, + "n_fits": len(fits), + "resolution_limited": resolution_limited, + } + + +def recommend_pid_gains(sysid: dict, guide_exp: float, verbose: bool, + lambda_l_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict: + """ + Recommend base RA/DEC proportional (+ conservative integral) gains from the + pulse_response step-response data, via a conservative SIMC/IMC-style PI rule + (pid_autotune_plan.md §3). Works for any mount type (WORM_GEAR, HARMONIC_DRIVE, + DIRECT_DRIVE) whose sysid data has pulse_response sessions; returns + confidence "unavailable" per axis when it doesn't. Calibration-time only, + never auto-applied -- see this module's docstring for how callers should + surface the result. + """ + return { + "ra": _recommend_axis_pid_gain(sysid, "RA", guide_exp, lambda_l_factor, verbose), + "dec": _recommend_axis_pid_gain(sysid, "DEC", guide_exp, lambda_l_factor, verbose), + } diff --git a/kstars/ekos/guide/offlinetrainer/pulse_response_fit.py b/kstars/ekos/guide/offlinetrainer/pulse_response_fit.py new file mode 100644 index 0000000000..c0dd9080f7 --- /dev/null +++ b/kstars/ekos/guide/offlinetrainer/pulse_response_fit.py @@ -0,0 +1,243 @@ +""" +offline_trainer/pulse_response_fit.py — generic step-response curve fitting for +pulse_response sysid sessions, shared across all three mount-type trainers. + +fit_pulse_response() has no mount-specific logic: it pairs opposite-direction +pulses (so PE/drift cancel in the difference), fits an exponential-approach +model per pulse magnitude/direction, and applies the same noise/sign- +consistency acceptance gates regardless of mount type. Two independent +consumers read its output: + + - train_harmonic.py uses (kappa, tau) directly as the Harmonic Drive + Kalman filter's spring/time-constant parameters. + - pid_autotune.py (all three mount types) uses the per-fit P_fit/tau_fit/ + residual_std/t_arr records (return_fits=True) as FOPDT step-response + data for the SIMC-style PID gain recommendation. + +Originally lived only in train_harmonic.py (the only mount type that +collected pulse_response data); extracted here once WORM_GEAR and +DIRECT_DRIVE gained their own pulse-response phases, since the fitting logic +was already 100% mount-agnostic. + +SPDX-License-Identifier: GPL-2.0-or-later +""" + +import numpy as np +import scipy.optimize + + +def fit_pulse_response(sysid: dict, axis: str, guide_exp: float, verbose: bool, + return_fits: bool = False): + """ + Fit spring constant κ and time constant τ from pulse_response sessions. + + Model: d(t) = P * (1 - κ * exp(-t/τ)) + v*t + c. Fits whose |P| is not + significantly above the residual noise are skipped. + + Returns: (kappa, tau_seconds), or (kappa, tau_seconds, fit_info) if + return_fits is True. fit_info["fits"] is the list of per-fit records + (P_fit_px, tau_fit_s, residual_std_px, t_first_s, t_arr, pos_arr) this + function already computes and would otherwise discard — reused by + pid_autotune.recommend_pid_gains() as the step-response data for PID + auto-tune, a second, independent consumer of the same pulse_response + sessions. fit_info["sign_consistent"] mirrors the gate this function + itself uses to decide the fits are real mechanics rather than noise. + """ + # Unmeasured means unmodeled: the default kappa stays 0 + DEFAULTS = (0.0, 1.5) + + def _finish(kappa, tau, fit_records, sign_consistent): + if not return_fits: + return kappa, tau + return kappa, tau, {"fits": fit_records, "sign_consistent": sign_consistent} + + pulse_sessions = [ + s for s in sysid["sessions"] + if s.get("type") == "pulse_response" and s.get("pulse_axis", "").upper() == axis.upper() + ] + + if not pulse_sessions: + if verbose: + print(f" [{axis}] No pulse_response sessions found. Using defaults (κ=0.2, τ=1.5s)") + return _finish(*DEFAULTS, [], None) + + axis_key = "ra_raw_px" if axis.upper() == "RA" else "dec_raw_px" + + def session_curve(s): + """(t, signed displacement from baseline) for one pulse session, or None.""" + frames = s.get("response_frames", []) + if len(frames) < 5: + return None + base = s.get("baseline_frames", []) + if base: + # New protocol: dedicated pre-pulse baseline; t is true seconds since the pulse. + baseline = float(np.mean([f.get(axis_key, 0.0) for f in base])) + t_vals = [f.get("t", (i + 1) * guide_exp) for i, f in enumerate(frames)] + pos_vals = [f.get(axis_key, 0.0) - baseline for f in frames] + else: + # Legacy: first frame doubles as the baseline. + baseline = frames[0].get(axis_key, 0.0) + t0 = frames[0].get("t", 0.0) + t_vals, pos_vals = [], [] + for i, f in enumerate(frames): + if i == 0: + continue + t = f.get("t", 0.0) - t0 + if t <= 0: + t = i * guide_exp + t_vals.append(t) + pos_vals.append(f.get(axis_key, 0.0) - baseline) + if len(t_vals) < 5: + return None + return np.array(t_vals, dtype=float), np.array(pos_vals, dtype=float) + + def fit_curve(t_arr, pos_arr, with_drift): + """Fit the spring model; returns (P, kappa, tau, residual_std) or None.""" + try: + if with_drift: + def model(t, P, kappa, tau, v, c): + return P * (1.0 - kappa * np.exp(-t / tau)) + v * t + c + slope0 = (pos_arr[-1] - pos_arr[0]) / max(t_arr[-1] - t_arr[0], 1e-3) + p0 = [pos_arr[-1] - slope0 * t_arr[-1], 0.3, 1.5, slope0, 0.0] + bounds = ([-50.0, 0.0, 0.1, -2.0, -10.0], [50.0, 0.9, 10.0, 2.0, 10.0]) + else: + def model(t, P, kappa, tau, c): + return P * (1.0 - kappa * np.exp(-t / tau)) + c + p0 = [pos_arr[-1], 0.3, 1.5, 0.0] + bounds = ([-100.0, 0.0, 0.1, -10.0], [100.0, 0.9, 10.0, 10.0]) + popt, _ = scipy.optimize.curve_fit(model, t_arr, pos_arr, p0=p0, + bounds=bounds, maxfev=10000) + residual_std = float(np.std(pos_arr - model(t_arr, *popt))) + return popt[0], popt[1], popt[2], residual_std + except (RuntimeError, ValueError): + return None + + kappas = [] + taus = [] + fit_signs = [] + paired_signs = set() + skipped_noise = 0 + fit_records = [] + + def accept_fit(kappa_fit, tau_fit, t_first): + # tau at the upper bound: exponential degenerate with the drift term + if tau_fit > 9.8: + return + # spring released before the first sample is indistinguishable from none + if tau_fit < t_first: + kappas.append(0.0) + else: + kappas.append(kappa_fit) + taus.append(tau_fit) + + # Pair opposite-direction sessions: the difference doubles the response + pos_dir, neg_dir = ("EAST", "WEST") if axis.upper() == "RA" else ("NORTH", "SOUTH") + by_mag = {} + for s in pulse_sessions: + by_mag.setdefault(s.get("pulse_magnitude_ms", 100.0), []).append(s) + + for pulse_mag, group in sorted(by_mag.items()): + pos_list = [s for s in group if s.get("pulse_direction", "").upper() == pos_dir] + neg_list = [s for s in group if s.get("pulse_direction", "").upper() == neg_dir] + paired = list(zip(pos_list, neg_list)) + leftovers = pos_list[len(paired):] + neg_list[len(paired):] + + for sp, sn in paired: + cp, cn = session_curve(sp), session_curve(sn) + if cp is None or cn is None: + continue + tp, pp = cp + tn, pn = cn + mask = (tp >= tn[0]) & (tp <= tn[-1]) + if mask.sum() < 5: + continue + t_arr = tp[mask] + diff = pp[mask] - np.interp(t_arr, tn, pn) + # Sessions are minutes apart so PE does not cancel exactly; v absorbs the leak + fit = fit_curve(t_arr, diff, with_drift=True) + if fit is None: + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms paired: curve_fit failed") + continue + P_fit, kappa_fit, tau_fit, residual_std = fit + if abs(P_fit) < 2.0 * residual_std: + skipped_noise += 1 + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms paired: |P|={abs(P_fit):.2f}px " + f"below noise ({residual_std:.2f}px) — skipped") + continue + paired_signs.add(1.0 if P_fit > 0 else -1.0) + accept_fit(kappa_fit, tau_fit, t_arr[0]) + fit_records.append({ + "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit), + "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std), + "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": diff, + }) + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms paired {pos_dir}-{neg_dir}: " + f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)") + + for s in leftovers: + c = session_curve(s) + if c is None: + continue + t_arr, pos_arr = c + fit = fit_curve(t_arr, pos_arr, with_drift=True) + if fit is None: + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms: curve_fit failed") + continue + P_fit, kappa_fit, tau_fit, residual_std = fit + if abs(P_fit) < 2.0 * residual_std: + skipped_noise += 1 + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms {s.get('pulse_direction', '?')}: " + f"response |P|={abs(P_fit):.2f}px below noise ({residual_std:.2f}px) — skipped") + continue + accept_fit(kappa_fit, tau_fit, t_arr[0]) + fit_signs.append((s.get("pulse_direction", "?"), np.sign(P_fit))) + fit_records.append({ + "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit), + "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std), + "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": pos_arr, + }) + if verbose: + print(f" [{axis}] Pulse {pulse_mag}ms {s.get('pulse_direction', '?')}: " + f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)") + + if not kappas: + if verbose: + print(f" [{axis}] No pulse response measurable above noise " + f"({skipped_noise} skipped). Using defaults (κ=0.2, τ=1.5s). " + f"Consider larger protocol pulses.") + return _finish(*DEFAULTS, fit_records, None) + + # Real responses have consistent signs per direction; paired diffs share one sign + by_dir = {} + for direction, sign in fit_signs: + by_dir.setdefault(direction, set()).add(sign) + dir_signs = [next(iter(s)) for s in by_dir.values() if len(s) == 1] + consistent = (len(paired_signs) <= 1 and + all(len(s) == 1 for s in by_dir.values()) and + (len(by_dir) < 2 or len(set(dir_signs)) == len(by_dir))) + if not consistent: + if verbose: + print(f" [{axis}] WARNING: response signs inconsistent across pulse directions " + f"— fits are noise, not mechanics. Using defaults (κ=0.2, τ=1.5s).") + return _finish(*DEFAULTS, fit_records, False) + + kappa_result = float(np.median(kappas)) + tau_result = float(np.median(taus)) if taus and kappa_result > 0.0 else DEFAULTS[1] + + # A median within ~2% of the fit bounds means the model chased noise/drift, not physics. + if kappa_result > 0.88 or tau_result > 9.8: + if verbose: + print(f" [{axis}] WARNING: fit pinned at bounds (κ={kappa_result:.3f}, " + f"τ={tau_result:.2f}s) — unphysical. Using defaults (κ=0.2, τ=1.5s).") + return _finish(*DEFAULTS, fit_records, consistent) + + if verbose: + print(f" [{axis}] Final: κ={kappa_result:.3f} (from {len(kappas)} fits), " + f"τ={tau_result:.2f}s") + + return _finish(kappa_result, tau_result, fit_records, consistent) diff --git a/kstars/ekos/guide/offlinetrainer/tail_ai_log.py b/kstars/ekos/guide/offlinetrainer/tail_ai_log.py new file mode 100644 index 0000000000..c317f4caf5 --- /dev/null +++ b/kstars/ekos/guide/offlinetrainer/tail_ai_log.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python3 +""" +tail_ai_log.py — Live tail/plot of the AI Guider debug CSV during a bench session. + +Watches the currently-growing ai_debug_logs/ai_guider_*.csv (KStars opens it once per +session and flushes a row every guide frame — see gmath.cpp's "AI DEBUG FILE LOGGER") +and plots a rolling window of: + 1. ra_error_arcsec / dec_error_arcsec — raw drift + 2. conf, with ai_state annotated on change + 3. ra_total_pulse_ms / dec_total_pulse_ms — the pulse actually sent to the mount, with + a zero line; this is the direct oscillation-watch panel + 4. ra_active_prop_gain — to correlate P-gain backoff with instability onset + +Panels 3-4 require the blend-breakdown columns added alongside this script; tailing an +older-format CSV degrades gracefully (those panels are skipped with a one-time warning). + +Usage: + python3 tail_ai_log.py # tails the newest file in the default dir + python3 tail_ai_log.py --log path/to/x.csv + python3 tail_ai_log.py --window 60 --refresh-ms 250 +""" + +import argparse +import csv +import glob +import os +import sys +from collections import deque + +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +DEFAULT_LOG_DIR = os.path.expanduser("~/.local/share/kstars/ai_debug_logs/") +# Hard cap so an all-night session doesn't grow memory unbounded; the plot only ever +# shows the last --window seconds anyway. +MAX_POINTS = 50000 + +BLEND_PROBE_COLUMNS = ("ra_total_pulse_ms", "dec_total_pulse_ms", "ra_active_prop_gain") + + +def find_latest_log(logdir: str) -> str: + candidates = glob.glob(os.path.join(logdir, "ai_guider_*.csv")) + if not candidates: + return None + return max(candidates, key=os.path.getmtime) + + +class LogTailer: + """Incrementally reads a growing CSV: tracks a byte offset and only parses newly + appended COMPLETE lines each poll, rather than re-reading the whole file.""" + + def __init__(self, path: str): + self.path = path + self.columns = None + self.col_index = {} + self.offset = 0 + self.buffers = {} + self.has_blend = False + self._warned_old_format = False + + def _init_header(self, f) -> bool: + header = f.readline() + if not header: + return False + self.columns = [c.strip() for c in header.decode("utf-8", errors="replace").strip().split(",")] + self.col_index = {name: i for i, name in enumerate(self.columns)} + self.has_blend = all(c in self.col_index for c in BLEND_PROBE_COLUMNS) + if not self.has_blend and not self._warned_old_format: + print("[warn] old-format CSV (no blend breakdown columns) — pulse/active-gain " + "panels will be blank.", file=sys.stderr) + self._warned_old_format = True + for name in self.columns: + self.buffers[name] = deque(maxlen=MAX_POINTS) + self.offset = f.tell() + return True + + def poll(self) -> bool: + """Read any newly-appended complete lines. Returns True if new rows were parsed.""" + if not os.path.exists(self.path): + return False + with open(self.path, "rb") as f: + if self.columns is None: + if not self._init_header(f): + return False + f.seek(self.offset) + raw = f.read() + if not raw: + return False + last_nl = raw.rfind(b"\n") + if last_nl < 0: + return False # partial line only; wait for the next poll + complete = raw[:last_nl + 1] + self.offset += len(complete) + + lines = complete.decode("utf-8", errors="replace").splitlines() + got_new = False + for row in csv.reader(lines): + if len(row) != len(self.columns): + continue + for name, val in zip(self.columns, row): + try: + self.buffers[name].append(float(val)) + except ValueError: + self.buffers[name].append(val) # non-numeric: algorithm/direction/state strings + got_new = True + return got_new + + def windowed(self, column: str, window_sec: float): + """Returns (t, values) for `column`, clipped to the last window_sec of t_session.""" + t = self.buffers.get("t_session") + v = self.buffers.get(column) + if not t or not v: + return [], [] + t_arr = list(t) + v_arr = list(v) + n = min(len(t_arr), len(v_arr)) + t_arr, v_arr = t_arr[-n:], v_arr[-n:] + if not t_arr: + return [], [] + t_max = t_arr[-1] + cutoff = t_max - window_sec + out_t, out_v = [], [] + for ti, vi in zip(t_arr, v_arr): + if ti >= cutoff: + out_t.append(ti) + out_v.append(vi) + return out_t, out_v + + +def main(): + parser = argparse.ArgumentParser(description="Live-tail and plot an AI Guider debug CSV.") + parser.add_argument("--log", type=str, default=None, + help="Path to the AI debug CSV (default: newest in " + f"{DEFAULT_LOG_DIR})") + parser.add_argument("--window", type=float, default=120.0, help="Rolling display window, seconds") + parser.add_argument("--refresh-ms", type=int, default=500, help="Poll/redraw interval, milliseconds") + args = parser.parse_args() + + log_path = args.log + if log_path is None: + log_path = find_latest_log(DEFAULT_LOG_DIR) + if log_path is None: + print(f"No log found yet in {DEFAULT_LOG_DIR} — waiting for guiding to start...") + + tailer = LogTailer(log_path) if log_path else None + + fig, axes = plt.subplots(4, 1, sharex=True, figsize=(11, 9)) + ax_err, ax_conf, ax_pulse, ax_gain = axes + fig.suptitle("AI Guider — live tail" + (f" ({os.path.basename(log_path)})" if log_path else " (waiting...)")) + + line_ra_err, = ax_err.plot([], [], label="RA error (arcsec)", color="#c0392b", linewidth=0.9) + line_dec_err, = ax_err.plot([], [], label="DEC error (arcsec)", color="#2980b9", linewidth=0.9) + ax_err.legend(loc="upper right", fontsize=8) + ax_err.set_ylabel("arcsec") + ax_err.grid(alpha=0.3) + + line_conf, = ax_conf.plot([], [], label="confidence", color="#27ae60", linewidth=1.2) + ax_conf.set_ylim(-0.05, 1.05) + ax_conf.set_ylabel("confidence") + ax_conf.grid(alpha=0.3) + ax_conf.legend(loc="upper right", fontsize=8) + + line_ra_pulse, = ax_pulse.plot([], [], label="RA total pulse (ms)", color="#c0392b", linewidth=0.9) + line_dec_pulse, = ax_pulse.plot([], [], label="DEC total pulse (ms)", color="#2980b9", linewidth=0.9) + ax_pulse.axhline(0, color="gray", linewidth=0.6) + ax_pulse.set_ylabel("ms") + ax_pulse.grid(alpha=0.3) + ax_pulse.legend(loc="upper right", fontsize=8) + + line_ra_gain, = ax_gain.plot([], [], label="RA active prop gain", color="#8e44ad", linewidth=0.9) + ax_gain.set_ylabel("gain") + ax_gain.set_xlabel("t_session (s)") + ax_gain.grid(alpha=0.3) + ax_gain.legend(loc="upper right", fontsize=8) + + state_annotations = {"last": None} + + def autoscale(ax, *y_lists): + vals = [v for lst in y_lists for v in lst if v == v] # drop NaN + if not vals: + return + lo, hi = min(vals), max(vals) + if lo == hi: + lo, hi = lo - 1, hi + 1 + pad = 0.1 * (hi - lo) + ax.set_ylim(lo - pad, hi + pad) + + def update(_frame): + nonlocal tailer, log_path + if tailer is None: + found = find_latest_log(DEFAULT_LOG_DIR) + if found is None: + return () + log_path = found + tailer = LogTailer(log_path) + fig.suptitle(f"AI Guider — live tail ({os.path.basename(log_path)})") + + tailer.poll() + w = args.window + + t_err, ra_err = tailer.windowed("ra_error_arcsec", w) + _, dec_err = tailer.windowed("dec_error_arcsec", w) + line_ra_err.set_data(t_err, ra_err) + line_dec_err.set_data(t_err, dec_err) + + t_conf, conf = tailer.windowed("conf", w) + line_conf.set_data(t_conf, conf) + + artists = [line_ra_err, line_dec_err, line_conf] + + if tailer.has_blend: + t_p, ra_pulse = tailer.windowed("ra_total_pulse_ms", w) + _, dec_pulse = tailer.windowed("dec_total_pulse_ms", w) + line_ra_pulse.set_data(t_p, ra_pulse) + line_dec_pulse.set_data(t_p, dec_pulse) + + t_g, ra_gain = tailer.windowed("ra_active_prop_gain", w) + line_ra_gain.set_data(t_g, ra_gain) + artists += [line_ra_pulse, line_dec_pulse, line_ra_gain] + + if t_err: + autoscale(ax_pulse, ra_pulse, dec_pulse) + autoscale(ax_gain, ra_gain) + + # Annotate ai_state transitions with a vertical marker. + ai_states = tailer.buffers.get("ai_state") + if ai_states and len(ai_states) > 0: + current_state = ai_states[-1] + if current_state != state_annotations["last"] and t_err: + ax_conf.axvline(t_err[-1], color="orange", linestyle="--", linewidth=0.8) + ax_conf.text(t_err[-1], 1.02, str(current_state), fontsize=7, color="orange", + ha="right", va="bottom", rotation=90) + state_annotations["last"] = current_state + + for ax in (ax_err, ax_conf, ax_pulse, ax_gain): + ax.relim() + ax.autoscale_view(scalex=True, scaley=False) + if t_err: + for ax in (ax_err, ax_conf, ax_pulse, ax_gain): + ax.set_xlim(max(0, t_err[0]), t_err[-1] + 0.5) + autoscale(ax_err, ra_err, dec_err) + + return artists + + anim = animation.FuncAnimation(fig, update, interval=args.refresh_ms, cache_frame_data=False) + plt.tight_layout() + plt.show() + del anim # keep a reference alive for the duration of plt.show(); silence unused warnings + + +if __name__ == "__main__": + main() diff --git a/kstars/ekos/guide/offlinetrainer/train.py b/kstars/ekos/guide/offlinetrainer/train.py index 9fa0662b9a..9fba8746c3 100644 --- a/kstars/ekos/guide/offlinetrainer/train.py +++ b/kstars/ekos/guide/offlinetrainer/train.py @@ -33,10 +33,10 @@ def parse_args(): p.add_argument("--epochs", type=int, default=None, help="Override default epoch count for neural models") p.add_argument("--pid-lambda-factor", type=float, default=None, - help="HARMONIC_DRIVE only: SIMC design parameter for the advisory PID " - "auto-tune recommendation (lambda = max(tau, factor * dead_time)). " - "Larger = slower/more conservative gain recommendation. " - "Default: train_harmonic.SIMC_LAMBDA_L_FACTOR") + help="SIMC design parameter for the advisory PID auto-tune recommendation " + "(lambda = max(tau, factor * dead_time)), computed for any mount type " + "whose sysid data includes pulse_response sessions. Larger = slower/" + "more conservative gain recommendation. Default: pid_autotune.SIMC_LAMBDA_L_FACTOR") p.add_argument("--simulate", action="store_true", help="Run closed-loop simulation after training completes") p.add_argument("--plot", action="store_true", @@ -77,20 +77,23 @@ def main(): print(f"[Ekos AI Trainer] Mount type: {mount_type}") print(f"[Ekos AI Trainer] Mount name: {sysid['equipment'].get('mount_name', 'unknown')}") + from pid_autotune import SIMC_LAMBDA_L_FACTOR + pid_lambda_factor = args.pid_lambda_factor or SIMC_LAMBDA_L_FACTOR + # Dispatch to the appropriate trainer if mount_type == "DIRECT_DRIVE": from train_direct_drive import train_direct_drive - weights = train_direct_drive(sysid, verbose=args.verbose) + weights = train_direct_drive(sysid, verbose=args.verbose, pid_lambda_factor=pid_lambda_factor) elif mount_type == "WORM_GEAR": from train_worm_gear import train_worm_gear - weights = train_worm_gear(sysid, gpu=args.gpu, - epochs=args.epochs, verbose=args.verbose) + weights = train_worm_gear(sysid, gpu=args.gpu, epochs=args.epochs, verbose=args.verbose, + pid_lambda_factor=pid_lambda_factor) elif mount_type == "HARMONIC_DRIVE": - from train_harmonic import train_harmonic, SIMC_LAMBDA_L_FACTOR + from train_harmonic import train_harmonic weights = train_harmonic(sysid, gpu=args.gpu, epochs=args.epochs, verbose=args.verbose, - pid_lambda_factor=args.pid_lambda_factor or SIMC_LAMBDA_L_FACTOR) + pid_lambda_factor=pid_lambda_factor) else: sys.exit( diff --git a/kstars/ekos/guide/offlinetrainer/train_direct_drive.py b/kstars/ekos/guide/offlinetrainer/train_direct_drive.py index 0d22dfcff6..c6b97c5ba8 100644 --- a/kstars/ekos/guide/offlinetrainer/train_direct_drive.py +++ b/kstars/ekos/guide/offlinetrainer/train_direct_drive.py @@ -18,6 +18,8 @@ import scipy.stats from typing import Optional from datetime import datetime +from pid_autotune import recommend_pid_gains, SIMC_LAMBDA_L_FACTOR + def _effective_pixel_scale(sysid): @@ -37,11 +39,20 @@ def _effective_pixel_scale(sysid): return ps def train_direct_drive(sysid: dict, - verbose: bool = False) -> dict: + verbose: bool = False, + pid_lambda_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict: """ Fit the 4-parameter refraction model from free-drift sysid sessions. Returns a weights dict compatible with DirectDriveGuider::loadWeights(). + + Also computes an advisory "pid_autotune" PID gain recommendation from any + pulse_response sessions present (see pid_autotune.py / pid_autotune_plan.md + §8). DIRECT_DRIVE mounts are expected to show a small, near-negligible tau/ + dead-time -- a confidently small result is itself a useful finding, not + just a null result. Returns confidence "unavailable" per axis if this + sysid run has no pulse_response data (Options::aIPIDAutoTune() was + off, or an older protocol run predates this mount type collecting it). """ eq = sysid["equipment"] pixel_scale = _effective_pixel_scale(sysid) # arcsec/px @@ -146,6 +157,14 @@ def train_direct_drive(sysid: dict, # Build model fingerprint from equipment block fingerprint = _build_fingerprint(sysid) + # Advisory PID auto-tune recommendation from pulse_response sessions, if any + # (see pid_autotune.py). Not applied automatically -- a human reviews and + # manually updates Options::rA/dECProportionalGain() in KStars. + pid_autotune = recommend_pid_gains(sysid, guide_exp, verbose, pid_lambda_factor) + + def _recommended(axis_result, field): + return axis_result[field] if axis_result["confidence"] != "unavailable" else None + return { "format_version": "1.0", "mount_type": "DIRECT_DRIVE", @@ -153,6 +172,13 @@ def train_direct_drive(sysid: dict, "mount_name": eq.get("mount_name", "unknown"), "pixel_scale": pixel_scale, "model_fingerprint": fingerprint, + # Advisory PID auto-tune recommendation -- see train_direct_drive()'s docstring. + # recommended_* fields are None when confidence is "unavailable" (don't apply). + "recommended_ra_proportional_gain": _recommended(pid_autotune["ra"], "proportional_gain"), + "recommended_ra_integral_gain": _recommended(pid_autotune["ra"], "integral_gain"), + "recommended_dec_proportional_gain": _recommended(pid_autotune["dec"], "proportional_gain"), + "recommended_dec_integral_gain": _recommended(pid_autotune["dec"], "integral_gain"), + "pid_autotune": pid_autotune, "parameters": { "k_ref": float(k_ref), "d_polar": float(d_polar), diff --git a/kstars/ekos/guide/offlinetrainer/train_harmonic.py b/kstars/ekos/guide/offlinetrainer/train_harmonic.py index 99a21d3017..5d946cc4b1 100644 --- a/kstars/ekos/guide/offlinetrainer/train_harmonic.py +++ b/kstars/ekos/guide/offlinetrainer/train_harmonic.py @@ -26,28 +26,10 @@ try: except ImportError: TORCH_AVAILABLE = False - - -# Default SIMC design parameter for the PID auto-tune recommendation (§ Phase 1b -# below): lambda = max(tau, SIMC_LAMBDA_L_FACTOR * dead_time). Larger -> slower/ -# more robust closed loop; smaller -> faster/less margin. Kept as an explicit, -# overridable constant per the plan's "expose lambda as a documented knob, not a -# hidden constant" (pid_autotune_plan.md §3.2). -SIMC_LAMBDA_L_FACTOR = 3.0 - -# KStars' "integral gain" (Options::rAIntegralGain()/dECIntegralGain()) multiplies -# a ~100-frame moving average of drift (gmath.cpp::processAxis(), drift_integral[k]), -# not a classical accumulating integrator. SIMC's tau_I (a reset *time*) has no -# principled mapping onto that EMA-style term, so rather than inventing one, the -# integral recommendation is a conservative fixed fraction of the proportional -# recommendation (pid_autotune_plan.md §1: "keep I conservative"). tau_I is still -# reported in the output for reference. -INTEGRAL_GAIN_CONSERVATIVE_FRACTION = 0.25 - -# Below this many usable step-response fits, the recommendation is flagged "low" -# confidence rather than withheld outright -- still informative, not to be applied -# unattended (pid_autotune_plan.md §4). -MIN_FITS_FOR_MEDIUM_CONFIDENCE = 6 +# Elastic-windup (kappa/tau) spring fitting is disabled for now -- see the note +# in train_harmonic() below. Re-import this if that fit is ever revisited. +# from pulse_response_fit import fit_pulse_response as _fit_spring_params +from pid_autotune import recommend_pid_gains as _recommend_pid_gains, SIMC_LAMBDA_L_FACTOR def _effective_pixel_scale(sysid): @@ -89,15 +71,29 @@ def train_harmonic(sysid: dict, guide_exp = eq.get("guide_exposure_ms", 1000.0) / 1000.0 if verbose: - print(f"\n--- Phase 1: Spring Parameter Fitting ---") - - # ── Step 1: Fit spring parameters from pulse_response sessions ───────── - kappa_ra, tau_ra = _fit_spring_params(sysid, "RA", guide_exp, verbose) - kappa_dec, tau_dec = _fit_spring_params(sysid, "DEC", guide_exp, verbose) + print(f"\n--- Phase 1: Spring Parameter Fitting (disabled -- see note below) ---") + + # ── Step 1: Elastic/spring (kappa/tau) fitting -- commented out for now ──── + # _fit_spring_params() (pulse_response_fit.fit_pulse_response()) has never + # resolved a spring response above the noise floor on any rig tested so far + # (pid_autotune_plan.md §9.1): kappa has come back 0.0 on every real run, + # bounded by the achieved guide-frame cadence being too coarse relative to + # any plausible spring time constant. Rather than spend the pulse_response + # phase's frame budget on a fit that has never once produced a usable + # result, this is disabled for now and left as a future-exploration item + # (e.g. revisit if a rig ever shows a resolvable spring, or add a separate + # opt-in toggle specifically for elastic-windup modeling rather than always + # attempting it as a side effect of PID auto-tune's data collection). + # The same pulse_response sessions still fully serve PID auto-tune (Step 1b + # below), which never depended on this fit succeeding. + kappa_ra, tau_ra = 0.0, 1.5 + kappa_dec, tau_dec = 0.0, 1.5 + # kappa_ra, tau_ra = _fit_spring_params(sysid, "RA", guide_exp, verbose) + # kappa_dec, tau_dec = _fit_spring_params(sysid, "DEC", guide_exp, verbose) if verbose: - print(f" κ_ra={kappa_ra:.3f} τ_ra={tau_ra:.2f}s") - print(f" κ_dec={kappa_dec:.3f} τ_dec={tau_dec:.2f}s") + print(f" κ_ra={kappa_ra:.3f} τ_ra={tau_ra:.2f}s (fit disabled, using default)") + print(f" κ_dec={kappa_dec:.3f} τ_dec={tau_dec:.2f}s (fit disabled, using default)") print(f"\n--- Phase 1b: PID Auto-tune (recommendation only, not applied) ---") # ── Step 1b: Recommend base P(+I) guide gains from the same step-response data ── @@ -174,355 +170,6 @@ def train_harmonic(sysid: dict, } -# ═══════════════════════════════════════════════════════════════════════════════ -# Phase 1: Spring parameter fitting from pulse_response sessions -# ═══════════════════════════════════════════════════════════════════════════════ - -def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool, - return_fits: bool = False): - """ - Fit spring constant κ and time constant τ from pulse_response sessions. - - Model: d(t) = P * (1 - κ * exp(-t/τ)) + v*t + c. Fits whose |P| is not - significantly above the residual noise are skipped. - - Returns: (kappa, tau_seconds), or (kappa, tau_seconds, fit_info) if - return_fits is True. fit_info["fits"] is the list of per-fit records - (P_fit_px, tau_fit_s, residual_std_px, t_first_s, t_arr, pos_arr) this - function already computes and would otherwise discard — reused by - _recommend_pid_gains() as the step-response data for PID auto-tune, a - second, independent consumer of the same pulse_response sessions. - fit_info["sign_consistent"] mirrors the gate this function itself uses - to decide the fits are real mechanics rather than noise. - """ - # Unmeasured means unmodeled: the default kappa stays 0 - DEFAULTS = (0.0, 1.5) - - def _finish(kappa, tau, fit_records, sign_consistent): - if not return_fits: - return kappa, tau - return kappa, tau, {"fits": fit_records, "sign_consistent": sign_consistent} - - pulse_sessions = [ - s for s in sysid["sessions"] - if s.get("type") == "pulse_response" and s.get("pulse_axis", "").upper() == axis.upper() - ] - - if not pulse_sessions: - if verbose: - print(f" [{axis}] No pulse_response sessions found. Using defaults (κ=0.2, τ=1.5s)") - return _finish(*DEFAULTS, [], None) - - axis_key = "ra_raw_px" if axis.upper() == "RA" else "dec_raw_px" - - def session_curve(s): - """(t, signed displacement from baseline) for one pulse session, or None.""" - frames = s.get("response_frames", []) - if len(frames) < 5: - return None - base = s.get("baseline_frames", []) - if base: - # New protocol: dedicated pre-pulse baseline; t is true seconds since the pulse. - baseline = float(np.mean([f.get(axis_key, 0.0) for f in base])) - t_vals = [f.get("t", (i + 1) * guide_exp) for i, f in enumerate(frames)] - pos_vals = [f.get(axis_key, 0.0) - baseline for f in frames] - else: - # Legacy: first frame doubles as the baseline. - baseline = frames[0].get(axis_key, 0.0) - t0 = frames[0].get("t", 0.0) - t_vals, pos_vals = [], [] - for i, f in enumerate(frames): - if i == 0: - continue - t = f.get("t", 0.0) - t0 - if t <= 0: - t = i * guide_exp - t_vals.append(t) - pos_vals.append(f.get(axis_key, 0.0) - baseline) - if len(t_vals) < 5: - return None - return np.array(t_vals, dtype=float), np.array(pos_vals, dtype=float) - - def fit_curve(t_arr, pos_arr, with_drift): - """Fit the spring model; returns (P, kappa, tau, residual_std) or None.""" - try: - if with_drift: - def model(t, P, kappa, tau, v, c): - return P * (1.0 - kappa * np.exp(-t / tau)) + v * t + c - slope0 = (pos_arr[-1] - pos_arr[0]) / max(t_arr[-1] - t_arr[0], 1e-3) - p0 = [pos_arr[-1] - slope0 * t_arr[-1], 0.3, 1.5, slope0, 0.0] - bounds = ([-50.0, 0.0, 0.1, -2.0, -10.0], [50.0, 0.9, 10.0, 2.0, 10.0]) - else: - def model(t, P, kappa, tau, c): - return P * (1.0 - kappa * np.exp(-t / tau)) + c - p0 = [pos_arr[-1], 0.3, 1.5, 0.0] - bounds = ([-100.0, 0.0, 0.1, -10.0], [100.0, 0.9, 10.0, 10.0]) - popt, _ = scipy.optimize.curve_fit(model, t_arr, pos_arr, p0=p0, - bounds=bounds, maxfev=10000) - residual_std = float(np.std(pos_arr - model(t_arr, *popt))) - return popt[0], popt[1], popt[2], residual_std - except (RuntimeError, ValueError): - return None - - kappas = [] - taus = [] - fit_signs = [] - paired_signs = set() - skipped_noise = 0 - fit_records = [] - - def accept_fit(kappa_fit, tau_fit, t_first): - # tau at the upper bound: exponential degenerate with the drift term - if tau_fit > 9.8: - return - # spring released before the first sample is indistinguishable from none - if tau_fit < t_first: - kappas.append(0.0) - else: - kappas.append(kappa_fit) - taus.append(tau_fit) - - # Pair opposite-direction sessions: the difference doubles the response - pos_dir, neg_dir = ("EAST", "WEST") if axis.upper() == "RA" else ("NORTH", "SOUTH") - by_mag = {} - for s in pulse_sessions: - by_mag.setdefault(s.get("pulse_magnitude_ms", 100.0), []).append(s) - - for pulse_mag, group in sorted(by_mag.items()): - pos_list = [s for s in group if s.get("pulse_direction", "").upper() == pos_dir] - neg_list = [s for s in group if s.get("pulse_direction", "").upper() == neg_dir] - paired = list(zip(pos_list, neg_list)) - leftovers = pos_list[len(paired):] + neg_list[len(paired):] - - for sp, sn in paired: - cp, cn = session_curve(sp), session_curve(sn) - if cp is None or cn is None: - continue - tp, pp = cp - tn, pn = cn - mask = (tp >= tn[0]) & (tp <= tn[-1]) - if mask.sum() < 5: - continue - t_arr = tp[mask] - diff = pp[mask] - np.interp(t_arr, tn, pn) - # Sessions are minutes apart so PE does not cancel exactly; v absorbs the leak - fit = fit_curve(t_arr, diff, with_drift=True) - if fit is None: - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms paired: curve_fit failed") - continue - P_fit, kappa_fit, tau_fit, residual_std = fit - if abs(P_fit) < 2.0 * residual_std: - skipped_noise += 1 - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms paired: |P|={abs(P_fit):.2f}px " - f"below noise ({residual_std:.2f}px) — skipped") - continue - paired_signs.add(1.0 if P_fit > 0 else -1.0) - accept_fit(kappa_fit, tau_fit, t_arr[0]) - fit_records.append({ - "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit), - "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std), - "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": diff, - }) - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms paired {pos_dir}-{neg_dir}: " - f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)") - - for s in leftovers: - c = session_curve(s) - if c is None: - continue - t_arr, pos_arr = c - fit = fit_curve(t_arr, pos_arr, with_drift=True) - if fit is None: - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms: curve_fit failed") - continue - P_fit, kappa_fit, tau_fit, residual_std = fit - if abs(P_fit) < 2.0 * residual_std: - skipped_noise += 1 - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms {s.get('pulse_direction', '?')}: " - f"response |P|={abs(P_fit):.2f}px below noise ({residual_std:.2f}px) — skipped") - continue - accept_fit(kappa_fit, tau_fit, t_arr[0]) - fit_signs.append((s.get("pulse_direction", "?"), np.sign(P_fit))) - fit_records.append({ - "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit), - "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std), - "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": pos_arr, - }) - if verbose: - print(f" [{axis}] Pulse {pulse_mag}ms {s.get('pulse_direction', '?')}: " - f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)") - - if not kappas: - if verbose: - print(f" [{axis}] No pulse response measurable above noise " - f"({skipped_noise} skipped). Using defaults (κ=0.2, τ=1.5s). " - f"Consider larger protocol pulses.") - return _finish(*DEFAULTS, fit_records, None) - - # Real responses have consistent signs per direction; paired diffs share one sign - by_dir = {} - for direction, sign in fit_signs: - by_dir.setdefault(direction, set()).add(sign) - dir_signs = [next(iter(s)) for s in by_dir.values() if len(s) == 1] - consistent = (len(paired_signs) <= 1 and - all(len(s) == 1 for s in by_dir.values()) and - (len(by_dir) < 2 or len(set(dir_signs)) == len(by_dir))) - if not consistent: - if verbose: - print(f" [{axis}] WARNING: response signs inconsistent across pulse directions " - f"— fits are noise, not mechanics. Using defaults (κ=0.2, τ=1.5s).") - return _finish(*DEFAULTS, fit_records, False) - - kappa_result = float(np.median(kappas)) - tau_result = float(np.median(taus)) if taus and kappa_result > 0.0 else DEFAULTS[1] - - # A median within ~2% of the fit bounds means the model chased noise/drift, not physics. - if kappa_result > 0.88 or tau_result > 9.8: - if verbose: - print(f" [{axis}] WARNING: fit pinned at bounds (κ={kappa_result:.3f}, " - f"τ={tau_result:.2f}s) — unphysical. Using defaults (κ=0.2, τ=1.5s).") - return _finish(*DEFAULTS, fit_records, consistent) - - if verbose: - print(f" [{axis}] Final: κ={kappa_result:.3f} (from {len(kappas)} fits), " - f"τ={tau_result:.2f}s") - - return _finish(kappa_result, tau_result, fit_records, consistent) - - -# ═══════════════════════════════════════════════════════════════════════════════ -# Phase 1b: PID auto-tune -- SIMC step-response gain recommendation -# -# Offline, calibration-time-only calculation (pid_autotune_plan.md). Reuses the -# per-fit P_fit/tau_fit/residual_std/t_arr values _fit_spring_params() already -# computes from pulse_response sessions -- no new data collection, just a second -# consumer of the existing step-response fits. -# ═══════════════════════════════════════════════════════════════════════════════ - -def _estimate_dead_time_s(t_arr: np.ndarray, pos_arr: np.ndarray, residual_std: float) -> float: - """ - First t at which |pos(t)| exceeds ~2.5x the fit's residual noise -- a proxy - for FOPDT dead time L (pid_autotune_plan.md §3.1). If even the first sample - already exceeds threshold (the common case at ~2-3s/frame cadence), this - returns t_arr[0]: L is only known to be <= the first sample, not resolved any - finer -- see _recommend_axis_pid_gain()'s "resolution_limited" flag. - """ - threshold = max(2.5 * residual_std, 1e-6) - for t, p in zip(t_arr, pos_arr): - if abs(p) > threshold: - return float(t) - return float(t_arr[-1]) - - -def _calibration_ms_per_arcsec(sysid: dict, axis: str) -> float: - """ - Median calibrated ms-per-arcsec pulse rate for this axis, from whichever - sessions recorded it (standard_guiding sessions carry ra_ms_per_arcsec/ - dec_ms_per_arcsec). This is Calibration::ra/decPulseMillisecondsPerArcsecond() - at collection time -- the same normalization gmath.cpp::processAxis() applies - to proportional_gain, needed to convert a physical SIMC Kc (ms/arcsec) back - into KStars' dimensionless aggressiveness (pid_autotune_plan.md §3.3). - """ - key = "ra_ms_per_arcsec" if axis.upper() == "RA" else "dec_ms_per_arcsec" - values = [s[key] for s in sysid["sessions"] if s.get(key, 0.0) and s[key] > 0.0] - return float(np.median(values)) if values else 0.0 - - -def _recommend_axis_pid_gain(sysid: dict, axis: str, guide_exp: float, - lambda_l_factor: float, verbose: bool) -> dict: - """ - Derive a recommended base proportional/integral gain for one axis from its - pulse_response step-response fits, via a conservative SIMC/IMC-style PI rule - (pid_autotune_plan.md §3.2). Returns a dict; confidence "unavailable" means - the numbers (if present at all) should not be used. - """ - _, _, fit_info = _fit_spring_params(sysid, axis, guide_exp, False, return_fits=True) - fits = fit_info["fits"] - - if fit_info["sign_consistent"] is False: - return {"confidence": "unavailable", - "reason": "pulse-response signs inconsistent across directions -- fits are noise, not mechanics"} - if not fits: - return {"confidence": "unavailable", "reason": "no usable pulse-response step-response fits"} - - cal_ms_per_arcsec = _calibration_ms_per_arcsec(sysid, axis) - if cal_ms_per_arcsec <= 0.0: - return {"confidence": "unavailable", - "reason": "no calibrated ms_per_arcsec recorded for this axis (need a standard_guiding session)"} - - pixel_scale = _effective_pixel_scale(sysid) - - K_samples, tau_samples, L_samples, t_first_samples = [], [], [], [] - for f in fits: - if f["pulse_magnitude_ms"] <= 0.0: - continue - K_samples.append(abs(f["P_fit_px"]) * pixel_scale / f["pulse_magnitude_ms"]) - t_first_samples.append(f["t_first_s"]) - L_samples.append(_estimate_dead_time_s(f["t_arr"], f["pos_arr"], f["residual_std_px"])) - if f["tau_fit_s"] <= 9.8: # same "pinned at bound, degenerate with drift" guard as the spring fit - tau_samples.append(f["tau_fit_s"]) - - if not K_samples: - return {"confidence": "unavailable", "reason": "no fit had a usable pulse magnitude"} - - K = float(np.median(K_samples)) # arcsec of steady-state response per ms of pulse - L = float(np.median(L_samples)) - # tau can't be resolved below the sampling floor either: a fit whose tau_fit - # landed below t_first ("spring already released") is only known to be - # <= t_first, not physically ~0 -- floor it the same way L is floored, so a - # spuriously tiny fitted tau can't produce a dangerously aggressive Kc. - tau_raw = float(np.median(tau_samples)) if tau_samples else float(np.median(t_first_samples)) - tau = max(tau_raw, L) - resolution_limited = bool(np.isclose(L, float(np.median(t_first_samples)), rtol=0.05)) - - lam = max(tau, lambda_l_factor * L) - Kc = (1.0 / K) * tau / (lam + L) # ms of pulse per arcsec of error - tau_i = min(tau, 4.0 * (lam + L)) # SIMC reset time, reported only -- see constant doc above - - proportional_gain = Kc / cal_ms_per_arcsec - integral_gain = INTEGRAL_GAIN_CONSERVATIVE_FRACTION * proportional_gain - confidence = "low" if (resolution_limited or len(fits) < MIN_FITS_FOR_MEDIUM_CONFIDENCE) else "medium" - - if verbose: - print(f" [{axis} PID] K={K:.5f} arcsec/ms tau={tau:.2f}s L={L:.2f}s lambda={lam:.2f}s " - f"(n={len(fits)} fits, cal={cal_ms_per_arcsec:.1f}ms/arcsec)") - print(f" [{axis} PID] Recommended proportional_gain={proportional_gain:.3f} " - f"integral_gain={integral_gain:.3f} confidence={confidence}") - - return { - "confidence": confidence, - "proportional_gain": float(np.clip(proportional_gain, 0.0, 1.0)), - "integral_gain": float(np.clip(integral_gain, 0.0, 1.0)), - "process_gain_arcsec_per_ms": K, - "tau_s": tau, - "dead_time_s": L, - "lambda_s": lam, - "tau_i_s": tau_i, - "calibration_ms_per_arcsec": cal_ms_per_arcsec, - "n_fits": len(fits), - "resolution_limited": resolution_limited, - } - - -def _recommend_pid_gains(sysid: dict, guide_exp: float, verbose: bool, - lambda_l_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict: - """ - Recommend base RA/DEC proportional (+ conservative integral) gains from the - pulse_response step-response data, via a conservative SIMC/IMC-style PI rule - (pid_autotune_plan.md §3). Calibration-time only, never auto-applied -- see - train_harmonic()'s docstring for how these surface in the weights JSON. - """ - return { - "ra": _recommend_axis_pid_gain(sysid, "RA", guide_exp, lambda_l_factor, verbose), - "dec": _recommend_axis_pid_gain(sysid, "DEC", guide_exp, lambda_l_factor, verbose), - } - - # ═══════════════════════════════════════════════════════════════════════════════ # Phase 2: PE period detection from free-drift data # ═══════════════════════════════════════════════════════════════════════════════ @@ -945,16 +592,17 @@ def _train_qnet(sysid, kappa_ra, tau_ra, kappa_dec, tau_dec, } -class QNet(nn.Module): - """5 → 8 (ReLU) → 2 Q-net for adaptive process noise.""" - def __init__(self): - super().__init__() - self.fc1 = nn.Linear(5, 8) - self.fc2 = nn.Linear(8, 2) +if TORCH_AVAILABLE: + class QNet(nn.Module): + """5 → 8 (ReLU) → 2 Q-net for adaptive process noise.""" + def __init__(self): + super().__init__() + self.fc1 = nn.Linear(5, 8) + self.fc2 = nn.Linear(8, 2) - def forward(self, x): - h = torch.relu(self.fc1(x)) - return self.fc2(h) + def forward(self, x): + h = torch.relu(self.fc1(x)) + return self.fc2(h) def _zero_qnet_weights(): diff --git a/kstars/ekos/guide/offlinetrainer/train_worm_gear.py b/kstars/ekos/guide/offlinetrainer/train_worm_gear.py index 293f83ca03..3c8b474d54 100644 --- a/kstars/ekos/guide/offlinetrainer/train_worm_gear.py +++ b/kstars/ekos/guide/offlinetrainer/train_worm_gear.py @@ -27,6 +27,8 @@ try: except ImportError: TORCH_AVAILABLE = False +from pid_autotune import recommend_pid_gains, SIMC_LAMBDA_L_FACTOR + def _effective_pixel_scale(sysid): @@ -48,10 +50,19 @@ def _effective_pixel_scale(sysid): def train_worm_gear(sysid: dict, gpu: bool = False, epochs: int = None, - verbose: bool = False) -> dict: + verbose: bool = False, + pid_lambda_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict: """ Train the PINN + residual MLP for a worm-gear mount. Returns a weights dict compatible with WormGearGuider::loadWeights(). + + Also computes an advisory "pid_autotune" PID gain recommendation from any + pulse_response sessions present (see pid_autotune.py / pid_autotune_plan.md + §8). Worm-gear mounts are exactly the class where backlash on DEC direction + reversal is a well-known effect -- it should show up directly as dead time + (L) in the FOPDT model. Returns confidence "unavailable" per axis if this + sysid run has no pulse_response data (Options::aIPIDAutoTune() was + off, or an older protocol run predates this mount type collecting it). """ if not TORCH_AVAILABLE: print("[ERROR] PyTorch is required to train the WormGearGuider MLP.") @@ -101,6 +112,14 @@ def train_worm_gear(sysid: dict, "dec_pulse_algorithm": fp.get("dec_pulse_algorithm", 0), } + # Advisory PID auto-tune recommendation from pulse_response sessions, if any + # (see pid_autotune.py). Not applied automatically -- a human reviews and + # manually updates Options::rA/dECProportionalGain() in KStars. + pid_autotune = recommend_pid_gains(sysid, guide_exp, verbose, pid_lambda_factor) + + def _recommended(axis_result, field): + return axis_result[field] if axis_result["confidence"] != "unavailable" else None + return { "format_version": "1.0", "mount_type": "WORM_GEAR", @@ -108,6 +127,13 @@ def train_worm_gear(sysid: dict, "mount_name": eq.get("mount_name", "unknown"), "pixel_scale": pixel_scale, "model_fingerprint": _build_fingerprint(sysid), + # Advisory PID auto-tune recommendation -- see train_worm_gear()'s docstring. + # recommended_* fields are None when confidence is "unavailable" (don't apply). + "recommended_ra_proportional_gain": _recommended(pid_autotune["ra"], "proportional_gain"), + "recommended_ra_integral_gain": _recommended(pid_autotune["ra"], "integral_gain"), + "recommended_dec_proportional_gain": _recommended(pid_autotune["dec"], "proportional_gain"), + "recommended_dec_integral_gain": _recommended(pid_autotune["dec"], "integral_gain"), + "pid_autotune": pid_autotune, "physics": { "pe_amplitude": float(pe_amplitude), "pe_period": float(pe_period), diff --git a/kstars/ekos/guide/offlinetrainer/validate_sysid.py b/kstars/ekos/guide/offlinetrainer/validate_sysid.py new file mode 100644 index 0000000000..72d1f679e4 --- /dev/null +++ b/kstars/ekos/guide/offlinetrainer/validate_sysid.py @@ -0,0 +1,336 @@ +#!/usr/bin/env python3 +""" +offline_trainer/validate_sysid.py — quick QC pass over a sysid_data_*.json file, +usable *during* a live AI Guiding Assistant / PID Auto-Tune run (the file is +rewritten after every completed session) or on a finished one, without waiting +for the full trainer. + +Written after a real bug (2026-08-03): once PID Auto-Tune was reordered to run +first in the protocol, pulse-response pulses could fire before guide calibration +had actually finished, silently recording ra_raw_px/dec_raw_px = 0.0 for every +single frame. Nothing about that looked wrong until someone loaded the JSON and +checked the raw numbers. This script exists to make that check a one-liner +instead of an ad hoc investigation. + +Usage: + python3 validate_sysid.py <sysid_data.json> [--verbose] + python3 validate_sysid.py --latest # newest file in ai_training_logs/ + python3 validate_sysid.py --latest --pid-lambda-factor 4.0 + +SPDX-License-Identifier: GPL-2.0-or-later +""" + +import argparse +import glob +import json +import os +import sys + +import numpy as np + +from pulse_response_fit import fit_pulse_response +from pid_autotune import recommend_pid_gains, SIMC_LAMBDA_L_FACTOR + +try: + from train_harmonic import _estimate_pe, _effective_pixel_scale + HAVE_PE = True +except ImportError: + HAVE_PE = False + + +AI_TRAINING_LOGS_DIR = os.path.expanduser("~/.local/share/kstars/ai_training_logs") + + +def _latest_sysid_file(): + files = sorted(glob.glob(os.path.join(AI_TRAINING_LOGS_DIR, "sysid_data_*.json"))) + if not files: + print(f"No sysid_data_*.json files found in {AI_TRAINING_LOGS_DIR}", file=sys.stderr) + sys.exit(1) + return files[-1] + + +def _pixel_scale(sysid): + if HAVE_PE: + return _effective_pixel_scale(sysid) + eq = sysid.get("equipment", {}) + return float(eq.get("pixel_scale_arcsec_per_px", 1.0) or 1.0) + + +def _rms_arcsec(frames, axis_key, pixel_scale): + vals = [f[axis_key] for f in frames if f.get("error_code", 0) == 0] + if not vals: + return None + return float(np.sqrt(np.mean(np.square(vals)))) * pixel_scale + + +def _check_pulse_response(sysid, guide_exp, pixel_scale, lambda_l_factor, verbose): + print("\n=== Pulse-response / PID Auto-Tune ===") + pr_sessions = [s for s in sysid["sessions"] if s.get("type") == "pulse_response"] + print(f" {len(pr_sessions)} pulse_response session(s) collected") + if not pr_sessions: + print(" (none yet -- nothing to check)") + return + + zero_sessions = 0 + for s in pr_sessions: + frames = s.get("response_frames", []) + if frames and all(f.get("ra_raw_px", 0.0) == 0.0 and f.get("dec_raw_px", 0.0) == 0.0 for f in frames): + zero_sessions += 1 + if zero_sessions: + print(f" *** WARNING: {zero_sessions}/{len(pr_sessions)} session(s) have ra_raw_px/dec_raw_px " + f"exactly 0.0 in EVERY frame. This is the calibration-not-ready bug (or the guider " + f"lost the star) -- pulses fired before real RA/DEC error was available. Don't trust " + f"any gain recommendation derived from this data.") + else: + print(" OK: no all-zero pulse_response sessions detected.") + + for axis in ("RA", "DEC"): + axis_sessions = [s for s in pr_sessions if s.get("pulse_axis") == axis] + if not axis_sessions: + continue + kappa, tau, info = fit_pulse_response(sysid, axis, guide_exp, verbose, return_fits=True) + fits = info["fits"] + print(f" [{axis}] {len(axis_sessions)} sessions -> {len(fits)} usable step-response fit(s), " + f"sign_consistent={info['sign_consistent']}") + if fits: + amps = [abs(f["P_fit_px"]) for f in fits] + print(f" amplitude range: {min(amps):.2f}-{max(amps):.2f} px " + f"(residual_std range: {min(f['residual_std_px'] for f in fits):.3f}" + f"-{max(f['residual_std_px'] for f in fits):.3f} px)") + + cal_note = ("(no calibrated ms/arcsec on record yet -- recommendation unavailable " + "until a standard_guiding session with ra/dec_ms_per_arcsec exists)") + gains = recommend_pid_gains(sysid, guide_exp, verbose, lambda_l_factor) + fp = sysid.get("model_fingerprint", {}) + for axis in ("ra", "dec"): + g = gains[axis] + if g["confidence"] == "unavailable": + print(f" [{axis.upper()}] recommendation unavailable: {g.get('reason', cal_note)}") + continue + res_flag = " *** resolution-limited (L pinned at first-sample time -- tau/L not really resolved)" \ + if g.get("resolution_limited") else "" + print(f" [{axis.upper()}] recommended proportional_gain={g['proportional_gain']:.3f} " + f"integral_gain={g['integral_gain']:.3f} (confidence={g['confidence']}, " + f"n_fits={g['n_fits']}, K={g['process_gain_arcsec_per_ms']:.5f}\"/ms, " + f"tau={g['tau_s']:.2f}s, L={g['dead_time_s']:.2f}s){res_flag}") + + # Cross-check: the live C++ gain-lock (a simplified plateau-average, no curve fit) + # may have already applied a gain -- compare it against this fuller offline recompute. + # A big disagreement is worth investigating (bad data, or the two implementations + # diverging) rather than assuming either one is automatically right. + live_gain = fp.get(f"{axis}_proportional_gain") + if live_gain is not None and g["proportional_gain"] > 0: + delta_pct = 100.0 * (live_gain - g["proportional_gain"]) / g["proportional_gain"] + if abs(delta_pct) > 15.0: + print(f" *** live-locked gain ({live_gain:.3f}) differs from this recompute " + f"by {delta_pct:+.0f}% -- worth checking why (C++ plateau-average vs. Python " + f"curve fit can differ some, but not usually this much)") + + +def _check_free_drift(sysid, pixel_scale, verbose): + print("\n=== Free-drift sessions ===") + sessions = [s for s in sysid["sessions"] if s.get("type") == "free_drift"] + print(f" {len(sessions)} free_drift session(s)") + for s in sessions: + frames = s.get("frames", []) + dur = s.get("duration_s", 0) + rms_ra = _rms_arcsec(frames, "ra_raw_px", pixel_scale) + rms_dec = _rms_arcsec(frames, "dec_raw_px", pixel_scale) + hf = s.get("hf_motion_arcsec") + flag = "" + if len(frames) < 10: + flag = " *** WARNING: very few frames -- segment likely too short to be useful" + print(f" {s.get('session_id', '?'):45s} dur={dur:5d}s frames={len(frames):4d} " + f"RMS(ra/dec)={rms_ra if rms_ra is None else round(rms_ra,2)}/" + f"{rms_dec if rms_dec is None else round(rms_dec,2)}\" " + f"hf_noise={hf if hf is None else round(hf,2)}\"{flag}") + + +def _check_standard_guiding(sysid, pixel_scale, verbose): + print("\n=== Standard-guiding sessions ===") + sessions = [s for s in sysid["sessions"] if s.get("type") == "standard_guiding"] + print(f" {len(sessions)} standard_guiding session(s)") + for s in sessions: + frames = s.get("frames", []) + dur = s.get("duration_s", 0) + rms_ra = _rms_arcsec(frames, "ra_raw_px", pixel_scale) + rms_dec = _rms_arcsec(frames, "dec_raw_px", pixel_scale) + gain_ra = s.get("aggressiveness_ra") + gain_dec = s.get("aggressiveness_dec") + lost = sum(1 for f in frames if f.get("error_code", 0) != 0) + flag = "" + if len(frames) < 10: + flag = " *** WARNING: very few frames -- segment likely too short to be useful" + elif lost > 0.1 * len(frames): + flag = f" *** WARNING: {lost}/{len(frames)} frames flagged star-lost" + print(f" {s.get('session_id', '?'):45s} dur={dur:5d}s frames={len(frames):4d} " + f"gain(ra/dec)={gain_ra}/{gain_dec} " + f"RMS(ra/dec)={rms_ra if rms_ra is None else round(rms_ra,2)}/" + f"{rms_dec if rms_dec is None else round(rms_dec,2)}\"{flag}") + + +def _check_artifacts(sysid, verbose): + """ + Cross-session anomaly/coverage summary -- the things worth flagging for improving + training data quality, the inference code, or the C++ AI guide engine itself, not + just per-session numbers. Returns a list of (severity, message) tuples for the + final verdict. + """ + print("\n=== Artifacts & coverage ===") + findings = [] + sessions = sysid.get("sessions", []) + + pr = [s for s in sessions if s.get("type") == "pulse_response"] + zero_pr = sum(1 for s in pr if s.get("response_frames") and + all(f.get("ra_raw_px", 0.0) == 0.0 and f.get("dec_raw_px", 0.0) == 0.0 + for f in s["response_frames"])) + if zero_pr: + findings.append(("HIGH", f"{zero_pr} pulse_response session(s) recorded all-zero " + f"raw_px -- calibration-not-ready or star-lost artifact")) + + # Duration shortfall: sessions that ran noticeably shorter than the phase asked for + # (dither/recenter interruptions, guider aborts) -- train.py silently accepts short + # segments, but they dilute the fit; worth knowing how much of the data is like this. + short_free_drift = short_std_guiding = 0 + for s in sessions: + if s.get("type") not in ("free_drift", "standard_guiding"): + continue + frames = s.get("frames", []) + # Rough expected frame count from duration_s / typical cadence isn't known here + # without guide_exp context per-call; just flag near-empty segments directly. + if len(frames) < 10: + if s["type"] == "free_drift": + short_free_drift += 1 + else: + short_std_guiding += 1 + if short_free_drift: + findings.append(("MEDIUM", f"{short_free_drift} free_drift segment(s) with <10 frames " + f"(interrupted early or barely started)")) + if short_std_guiding: + findings.append(("MEDIUM", f"{short_std_guiding} standard_guiding segment(s) with <10 frames")) + + # Star-lost rate across all non-pulse sessions. + total_frames = total_lost = 0 + for s in sessions: + if s.get("type") not in ("free_drift", "standard_guiding"): + continue + for f in s.get("frames", []): + total_frames += 1 + if f.get("error_code", 0) != 0: + total_lost += 1 + if total_frames > 0: + lost_pct = 100.0 * total_lost / total_frames + print(f" Star-lost frames: {total_lost}/{total_frames} ({lost_pct:.1f}%)") + if lost_pct > 5.0: + findings.append(("MEDIUM", f"{lost_pct:.1f}% of frames flagged star-lost across the run " + f"-- check guide star SNR/exposure for this rig")) + + # Pier-side and sky coverage -- directly relevant both to training (drift/refraction fit + # validity range) and to planning the later live-evaluation position matrix. + pier_sides = {s.get("pier_side") for s in sessions if s.get("pier_side")} + alts = [s.get("altitude_deg") for s in sessions if s.get("altitude_deg") is not None] + azs = [s.get("azimuth_deg") for s in sessions if s.get("azimuth_deg") is not None] + print(f" Pier sides tested: {sorted(pier_sides) if pier_sides else 'unknown'}") + if alts: + print(f" Altitude range covered: {min(alts):.0f}-{max(alts):.0f} deg") + if len(pier_sides) < 2: + findings.append(("INFO", f"Only pier side(s) {sorted(pier_sides)} tested in this data collection " + f"run -- the live evaluation phase should deliberately cover the other " + f"side (meridian flip) since drift/refraction/backlash can differ.")) + + # PE harmonic count vs. the live HarmonicGuider's fixed 2-line cap. + pe_lines = sysid.get("physical", {}).get("pe_lines") if "physical" in sysid else None + # sysid files from a live (not-yet-trained) run don't have "physical" -- that's only in + # a trained weights.json. Nothing to check here pre-training; see the trained-weights + # inspection step in the test plan instead. + + noise_floor = sysid.get("noise_floor_arcsec") + if noise_floor is not None: + print(f" Recorded seeing/noise floor (unguided HF star motion): {noise_floor:.2f}\"") + if noise_floor > 2.0: + findings.append(("INFO", f"Seeing/noise floor is {noise_floor:.2f}\" -- unusually high; " + f"any RMS close to this floor is at the limit of what guiding " + f"can physically improve, not a guider defect")) + else: + findings.append(("INFO", "No noise_floor_arcsec recorded yet (needs a free_drift session " + ">30 frames to measure)")) + + if not findings: + print(" No artifacts detected.") + else: + for severity, msg in findings: + print(f" [{severity}] {msg}") + return findings + + +def _check_pe(sysid, guide_exp, verbose): + if not HAVE_PE: + return + print("\n=== PE period preview (train_harmonic._estimate_pe) ===") + period, amplitude, lines = _estimate_pe(sysid, guide_exp, verbose) + if period <= 0: + print(" No significant PE detected yet (may need more free_drift/standard_guiding data).") + return + print(f" Primary period: {period:.1f}s amplitude: {amplitude:.3f}px") + for line in lines[1:]: + print(f" Secondary: {line['period_s']:.1f}s amplitude={line['amplitude_px']:.3f}px SNR={line['snr']:.0f}") + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("sysid_file", nargs="?", help="Path to sysid_data_*.json") + ap.add_argument("--latest", action="store_true", help="Use the newest file in ai_training_logs/") + ap.add_argument("--pid-lambda-factor", type=float, default=SIMC_LAMBDA_L_FACTOR) + ap.add_argument("--verbose", action="store_true") + args = ap.parse_args() + + if args.latest or not args.sysid_file: + path = _latest_sysid_file() + else: + path = args.sysid_file + + print(f"Loading {path}") + with open(path) as f: + sysid = json.load(f) + + eq = sysid.get("equipment", {}) + fp = sysid.get("model_fingerprint", {}) + guide_exp = eq.get("guide_exposure_ms", 1000.0) / 1000.0 + pixel_scale = _pixel_scale(sysid) + + print(f"\nmount_type={eq.get('mount_type')} mount_name={eq.get('mount_name')} " + f"camera={eq.get('camera')} pixel_scale={pixel_scale:.4f}\"/px guide_exp={guide_exp:.2f}s") + print(f"Current fingerprint: ra_gain={fp.get('ra_proportional_gain')} " + f"dec_gain={fp.get('dec_proportional_gain')} " + f"ra_int={fp.get('ra_integral_gain')} dec_int={fp.get('dec_integral_gain')}") + + from collections import Counter + counts = Counter(s.get("type") for s in sysid.get("sessions", [])) + print(f"Sessions so far: {dict(counts)}") + + _check_pulse_response(sysid, guide_exp, pixel_scale, args.pid_lambda_factor, args.verbose) + _check_free_drift(sysid, pixel_scale, args.verbose) + _check_standard_guiding(sysid, pixel_scale, args.verbose) + findings = _check_artifacts(sysid, args.verbose) + _check_pe(sysid, guide_exp, args.verbose) + + print("\n=== Verdict ===") + high = [m for sev, m in findings if sev == "HIGH"] + medium = [m for sev, m in findings if sev == "MEDIUM"] + if high: + print(" NOT READY -- fix the HIGH-severity issue(s) above before training or trusting any " + "gain recommendation from this data:") + for m in high: + print(f" - {m}") + elif medium: + print(" USABLE WITH CAVEATS -- training will run, but review the MEDIUM findings above:") + for m in medium: + print(f" - {m}") + else: + print(" Looks OK -- no high/medium-severity issues detected.") + + print("\nDone.") + + +if __name__ == "__main__": + main()