[education/kstars] kstars/ekos/guide: AI Guide: guard the drift fit against short segments, clamp parallactic extrapolation, Add debug logs
Jasem Mutlaq <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit 989665604dd2346d027e77593298f7185e581b82 by Jasem Mutlaq, on behalf of Pavan Kumar S G.
Committed on 05/08/2026 at 14:02.
Pushed by mutlaqja into branch 'master'.
AI Guide: guard the drift fit against short segments, clamp parallactic extrapolation, Add debug logs
This reverts commit f1f5ee41f8cd98a844f21ea392fb3da44f9314f2.
guard the drift fit against short segments, clamp parallactic extrapolation, Add debug logs
M +7 -1 kstars/ekos/guide/aiguideprotocol.cpp
M +13 -7 kstars/ekos/guide/internalguide/gmath.cpp
M +32 -21 kstars/ekos/guide/internalguide/harmonic_guider.cpp
M +2 -4 kstars/ekos/guide/internalguide/harmonic_guider.h
M +2 -0 kstars/ekos/guide/internalguide/mount_guider.h
M +35 -14 kstars/ekos/guide/offlinetrainer/train_harmonic.py
https://invent.kde.org/education/kstars/-/commit/989665604dd2346d027e77593298f7185e581b82
diff --git a/kstars/ekos/guide/aiguideprotocol.cpp b/kstars/ekos/guide/aiguideprotocol.cpp
index 50c89e3422..5678ca2a42 100644
--- a/kstars/ekos/guide/aiguideprotocol.cpp
+++ b/kstars/ekos/guide/aiguideprotocol.cpp
@@ -818,7 +818,13 @@ void AIGuideProtocol::processProtocol()
// trainer's truncation guard can still discard them.
const int recordedDuration = (phase.freeDrift && !m_PhaseAborted)
? m_SegmentSeconds : phase.durationSeconds;
- flushPhaseSegment(phase, recordedDuration);
+ // A drift tail too short to fit is not worth recording — the same rule the
+ // re-center path applies (a 12-frame tail once set k_ref_dec from noise).
+ if (phase.freeDrift && !m_PhaseAborted && m_PhaseData.size() < MIN_SEGMENT_FRAMES)
+ emit protocolLog(QString("Final drift segment too short (%1 frames) — discarded.")
+ .arg(m_PhaseData.size()));
+ else
+ flushPhaseSegment(phase, recordedDuration);
m_Phases.removeFirst();
m_State = STATE_PRECHECK;
diff --git a/kstars/ekos/guide/internalguide/gmath.cpp b/kstars/ekos/guide/internalguide/gmath.cpp
index ae30e4882a..c442457cc2 100644
--- a/kstars/ekos/guide/internalguide/gmath.cpp
+++ b/kstars/ekos/guide/internalguide/gmath.cpp
@@ -277,7 +277,7 @@ void cgmath::start()
{
qCInfo(KSTARS_EKOS_GUIDE) << "[AI GUIDER] Applied recorded settings:" << m_AIGuider->fingerprintApplied();
emit newLog(i18n("AI Guider: applied the trained model's recorded settings (%1).",
- m_AIGuider->fingerprintApplied()));
+ m_AIGuider->fingerprintApplied()));
}
if (useAIAlgorithm)
{
@@ -369,14 +369,14 @@ bool cgmath::reloadAIWeights()
const QString reason = candidate ? candidate->fingerprintError() : QString();
qCWarning(KSTARS_EKOS_GUIDE) << ">>> AI GUIDER RELOAD FAILED:" << weightsPath << reason;
emit newLog(i18n("AI Guider weight reload failed: %1", reason.isEmpty()
- ? i18n("the weights file could not be read or does not match this mount type.") : reason));
+ ? i18n("the weights file could not be read or does not match this mount type.") : reason));
return false;
}
m_AIGuider = std::move(candidate);
if (!m_AIGuider->fingerprintApplied().isEmpty())
emit newLog(i18n("AI Guider: applied the trained model's recorded settings (%1).",
- m_AIGuider->fingerprintApplied()));
+ m_AIGuider->fingerprintApplied()));
// A plain (non-forced) reset is enough: HarmonicGuider::resetSession() already forces a
// full reset itself when the new weights describe a different PE period than the currently
// active static Kalman state, and otherwise only clears transient tracking state while
@@ -1267,12 +1267,16 @@ void cgmath::performProcessing(Ekos::GuideState state, QSharedPointer<FITSData>
{
out << "t_session,dt,altitude_deg,azimuth_deg,parallactic_angle_deg,ra_error_arcsec,uncorrected_ra_delta_px,dec_error_arcsec,uncorrected_dec_delta_px,conf,pred_ra_arcsec,physics_ra_arcsec,mlp_ra_arcsec,pred_dec_arcsec,physics_dec_arcsec,mlp_dec_arcsec,ai_state,pe_statestring,"
"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\n";
+ "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,"
+ "ra_drift_arcsec,dec_drift_arcsec\n";
m_AIDebugHeaderWritten = true;
}
// Empty field for NaN (fields not meaningful for the algorithm that actually ran)
// rather than the literal string "nan", which not every CSV reader parses cleanly.
- auto csvNum = [](double v) { return std::isnan(v) ? QString() : QString::number(v, 'f', 3); };
+ auto csvNum = [](double v)
+ {
+ return std::isnan(v) ? QString() : QString::number(v, 'f', 3);
+ };
// Signed total pulse actually sent this frame, mirroring the sign convention
// updateOutParams() uses for m_accumulated_pulse_ra/dec (gmath.cpp:signed_pulse).
auto axisTotalMs = [this](int k) -> double
@@ -1280,7 +1284,7 @@ void cgmath::performProcessing(Ekos::GuideState state, QSharedPointer<FITSData>
if (out_params.pulse_dir[k] == NO_DIR || out_params.pulse_length[k] == 0)
return 0.0;
const bool neg = (k == GUIDE_RA && out_params.pulse_dir[k] == RA_DEC_DIR)
- || (k == GUIDE_DEC && out_params.pulse_dir[k] == DEC_DEC_DIR);
+ || (k == GUIDE_DEC && out_params.pulse_dir[k] == DEC_DEC_DIR);
return neg ? -out_params.pulse_length[k] : out_params.pulse_length[k];
};
out << frameData.t_session_sec << ","
@@ -1316,7 +1320,9 @@ void cgmath::performProcessing(Ekos::GuideState state, QSharedPointer<FITSData>
<< csvNum(m_lastBlend[GUIDE_DEC].activePropGain) << ","
<< axisTotalMs(GUIDE_DEC) << ","
<< directionStr(out_params.pulse_dir[GUIDE_DEC]) << ","
- << (out_params.pulse_dir[GUIDE_DEC] == NO_DIR ? 1 : 0) << "\n";
+ << (out_params.pulse_dir[GUIDE_DEC] == NO_DIR ? 1 : 0) << ","
+ << m_lastAIPrediction.drift_ra_arcsec << ","
+ << m_lastAIPrediction.drift_dec_arcsec << "\n";
out.flush();
}
}
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.cpp b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
index e11549dc31..d80e783ac2 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.cpp
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
@@ -114,6 +114,8 @@ bool HarmonicGuider::loadWeights(const QString &weightsPath)
}
m_fit_alt_min = phys["fit_alt_min"].toDouble(35.0);
m_fit_alt_max = phys["fit_alt_max"].toDouble(65.0);
+ m_fit_par_min = phys["fit_par_min"].toDouble(-90.0);
+ m_fit_par_max = phys["fit_par_max"].toDouble(90.0);
// Sanity bounds on spring parameters
m_kappa_ra = std::clamp(m_kappa_ra, 0.0, 0.9);
@@ -313,22 +315,20 @@ HarmonicGuider::computeQ(double snr, double snr_delta,
Q(RA_SPRING, RA_SPRING) = 0.001 * dt;
Q(DEC_SPRING, DEC_SPRING) = 0.001 * dt;
- // PE process noise (small — PE is nearly deterministic). RA only: strain-wave
- // periodic error comes from the continuously rotating RA drive, so a DEC PE
- // oscillator has nothing real to lock onto. With no measurement on those states
- // either (see kalmanUpdate()), a slow sinusoid and DEC_VEL's linear ramp become
- // degenerate over a session shorter than one PE period, and the filter settles
- // into a large, growing, mutually-cancelling pair whose small residual is a
- // steadily growing DEC bias.
+ // PE process noise (small — PE is nearly deterministic)
if (m_pe_period > 0.0)
{
Q(RA_PE_SIN, RA_PE_SIN) = 0.001 * dt;
Q(RA_PE_COS, RA_PE_COS) = 0.001 * dt;
+ Q(DEC_PE_SIN, DEC_PE_SIN) = 0.001 * dt;
+ Q(DEC_PE_COS, DEC_PE_COS) = 0.001 * dt;
}
if (m_pe2_period > 0.0)
{
Q(RA_PE2_SIN, RA_PE2_SIN) = 0.001 * dt;
Q(RA_PE2_COS, RA_PE2_COS) = 0.001 * dt;
+ Q(DEC_PE2_SIN, DEC_PE2_SIN) = 0.001 * dt;
+ Q(DEC_PE2_COS, DEC_PE2_COS) = 0.001 * dt;
}
return Q;
@@ -340,7 +340,10 @@ void HarmonicGuider::driftRates(double alt_deg, double parallactic_angle_deg,
{
const double alt_rad = std::clamp(alt_deg, m_fit_alt_min, m_fit_alt_max) * M_PI / 180.0;
const double cos_alt = std::cos(alt_rad);
- const double q_rad = parallactic_angle_deg * M_PI / 180.0;
+ // Clamp to the fitted range like altitude: a coefficient fitted at sin(q) <= 0.1 must
+ // not be extrapolated to sin(q) = 1 (that turned a noise fit into a -0.32"/s phantom
+ // DEC drift on real data — see DEC_OFFSET_ROOT_CAUSE.md).
+ const double q_rad = std::clamp(parallactic_angle_deg, m_fit_par_min, m_fit_par_max) * M_PI / 180.0;
ra_rate = m_drift_ra;
dec_rate = m_drift_dec + m_d_polar;
@@ -387,20 +390,20 @@ void HarmonicGuider::kalmanPredict(double dt, double alt_deg, double parallactic
// ── Kalman update step ───────────────────────────────────────────────────────
void HarmonicGuider::kalmanUpdate(double ra_meas_px, double dec_meas_px, double snr)
{
- // Observation matrix: observe position + PE_sin, RA only (see computeQ() for why
- // DEC has no PE oscillator to observe). DEC_PE_SIN/DEC_PE2_SIN never receive a
- // Kalman gain and stay at their zero-initialized value for the whole session.
- // H extracts: ra_obs = ra_err + pe_sin_ra, dec_obs = dec_err
+ // Observation matrix: observe position + PE_sin
+ // H extracts: ra_obs = ra_err + pe_sin_ra, dec_obs = dec_err + pe_sin_dec
Eigen::Matrix<double, N_OBS, N_STATES> H = Eigen::Matrix<double, N_OBS, N_STATES>::Zero();
H(0, RA_POS) = 1.0;
H(1, DEC_POS) = 1.0;
if (m_pe_period > 0.0)
{
H(0, RA_PE_SIN) = 1.0;
+ H(1, DEC_PE_SIN) = 1.0;
}
if (m_pe2_period > 0.0)
{
H(0, RA_PE2_SIN) = 1.0;
+ H(1, DEC_PE2_SIN) = 1.0;
}
// Measurement noise from the current frame's SNR (~0.5 px at SNR 30)
@@ -512,13 +515,19 @@ GuideOutput HarmonicGuider::predict(const GuideFrameData &frame)
out.ra_correction_arcsec = m_lastPredRA * frame.pixel_scale;
out.dec_correction_arcsec = m_lastPredDEC * frame.pixel_scale;
- // Debug breakdown: physics = trend (VEL), mlp = spring + PE (the "learned" part)
- const double drift_ra_px = m_x(RA_VEL) * frame.dt;
- const double drift_dec_px = m_x(DEC_VEL) * frame.dt;
- out.physics_ra_arcsec = drift_ra_px * frame.pixel_scale;
- out.physics_dec_arcsec = drift_dec_px * frame.pixel_scale;
- out.mlp_ra_arcsec = (m_lastPredRA - drift_ra_px) * frame.pixel_scale;
- out.mlp_dec_arcsec = (m_lastPredDEC - drift_dec_px) * frame.pixel_scale;
+ // Debug breakdown: physics = trend (VEL), drift = static model injection,
+ // mlp = spring + PE. Keeping the injection out of mlp_* matters: bundled together
+ // they made a bad drift coefficient look like an oscillator fault.
+ const double vel_ra_px = m_x(RA_VEL) * frame.dt;
+ const double vel_dec_px = m_x(DEC_VEL) * frame.dt;
+ double inj_ra_rate = 0.0, inj_dec_rate = 0.0;
+ driftRates(frame.altitude_deg, frame.parallactic_angle_deg, inj_ra_rate, inj_dec_rate);
+ out.physics_ra_arcsec = vel_ra_px * frame.pixel_scale;
+ out.physics_dec_arcsec = vel_dec_px * frame.pixel_scale;
+ out.drift_ra_arcsec = inj_ra_rate * frame.dt * frame.pixel_scale;
+ out.drift_dec_arcsec = inj_dec_rate * frame.dt * frame.pixel_scale;
+ out.mlp_ra_arcsec = (m_lastPredRA - vel_ra_px - inj_ra_rate * frame.dt) * frame.pixel_scale;
+ out.mlp_dec_arcsec = (m_lastPredDEC - vel_dec_px - inj_dec_rate * frame.dt) * frame.pixel_scale;
return out;
}
@@ -660,8 +669,10 @@ GuideOutput HarmonicGuider::darkPredict(double dt_sec)
out.physics_ra_arcsec = x(RA_VEL) * dt_sec * m_lastPixelScale;
out.physics_dec_arcsec = x(DEC_VEL) * dt_sec * m_lastPixelScale;
- out.mlp_ra_arcsec = (pred_ra - post_ra - x(RA_VEL) * dt_sec) * m_lastPixelScale;
- out.mlp_dec_arcsec = (pred_dec - post_dec - x(DEC_VEL) * dt_sec) * m_lastPixelScale;
+ out.drift_ra_arcsec = ra_rate * dt_sec * m_lastPixelScale;
+ out.drift_dec_arcsec = dec_rate * dt_sec * m_lastPixelScale;
+ out.mlp_ra_arcsec = (pred_ra - post_ra - (x(RA_VEL) + ra_rate) * dt_sec) * m_lastPixelScale;
+ out.mlp_dec_arcsec = (pred_dec - post_dec - (x(DEC_VEL) + dec_rate) * dt_sec) * m_lastPixelScale;
return out;
}
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.h b/kstars/ekos/guide/internalguide/harmonic_guider.h
index 4bee0ab68a..a7b3f71583 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.h
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.h
@@ -15,10 +15,6 @@
*
* The PE states evolve as a 2D rotation at the detected PE frequency, allowing
* the Kalman filter to automatically estimate PE amplitude and phase online.
- * PE is an RA-only phenomenon (strain-wave error from the continuously rotating
- * RA drive); pe_sin_dec/pe_cos_dec/pe2_sin_dec/pe2_cos_dec occupy state-vector
- * slots but are never driven by process noise or observed (see computeQ() and
- * kalmanUpdate()), so they stay at zero for the whole session.
*
* SPDX-License-Identifier: GPL-2.0-or-later
*/
@@ -96,6 +92,8 @@ class HarmonicGuider : public MountSpecificGuider
double m_k_ref_dec { 0.0 }; ///< DEC refraction coefficient
double m_fit_alt_min { 35.0 }; ///< Altitude range the drift fit is valid for
double m_fit_alt_max { 65.0 };
+ double m_fit_par_min { -90.0 }; ///< Parallactic range the DEC drift fit is valid for
+ double m_fit_par_max { 90.0 };
// ── Kalman filter state (static to survive object recreation) ────────
static Eigen::Matrix<double, N_STATES, 1> m_x; ///< State estimate
diff --git a/kstars/ekos/guide/internalguide/mount_guider.h b/kstars/ekos/guide/internalguide/mount_guider.h
index 4850fbdbc8..f633412e9c 100644
--- a/kstars/ekos/guide/internalguide/mount_guider.h
+++ b/kstars/ekos/guide/internalguide/mount_guider.h
@@ -67,6 +67,8 @@ struct GuideOutput
double dec_correction_arcsec { 0.0 }; ///< Feed-forward DEC correction (arcsec)
// Debug specific outputs
+ double drift_ra_arcsec { 0.0 }; ///< Static drift-model injection (RA)
+ double drift_dec_arcsec { 0.0 }; ///< Static drift-model injection (DEC)
double physics_ra_arcsec { 0.0 }; ///< Physics layer RA drift prediction
double physics_dec_arcsec { 0.0 }; ///< Physics layer DEC drift prediction
double mlp_ra_arcsec { 0.0 }; ///< MLP layer RA residual prediction
diff --git a/kstars/ekos/guide/offlinetrainer/train_harmonic.py b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
index 5d946cc4b1..f3970e3e10 100644
--- a/kstars/ekos/guide/offlinetrainer/train_harmonic.py
+++ b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
@@ -113,14 +113,15 @@ def train_harmonic(sysid: dict,
print(f"\n--- Phase 3: Drift Parameter Fitting ---")
# ── Step 3: Fit drift parameters from free_drift sessions ──────────────
- drift_ra, drift_dec, d_polar, k_ref, k_ref_dec = _fit_drift_params(
- sysid, guide_exp, verbose)
+ (drift_ra, drift_dec, d_polar, k_ref, k_ref_dec,
+ fit_alts, fit_qs) = _fit_drift_params(sysid, guide_exp, verbose)
- # Altitude range the drift/refraction fit is valid for (runtime clamps to it).
- fit_alts = [s.get("altitude_deg", 45.0) for s in sysid["sessions"]
- if s.get("type") == "free_drift" and len(s.get("frames", [])) >= 10]
+ # Geometry the drift/refraction fit is valid for (runtime clamps to it) — from the
+ # sessions the fit actually used, not everything recorded.
fit_alt_min = min(fit_alts) if fit_alts else 35.0
fit_alt_max = max(fit_alts) if fit_alts else 65.0
+ fit_par_min = min(fit_qs) if fit_qs else -90.0
+ fit_par_max = max(fit_qs) if fit_qs else 90.0
if verbose:
print(f" drift_ra={drift_ra:.6e} px/s drift_dec={drift_dec:.6e} px/s")
@@ -165,6 +166,8 @@ def train_harmonic(sysid: dict,
"k_ref_dec": float(k_ref_dec),
"fit_alt_min": float(fit_alt_min),
"fit_alt_max": float(fit_alt_max),
+ "fit_par_min": float(fit_par_min),
+ "fit_par_max": float(fit_par_max),
},
"qnet": qnet_weights,
}
@@ -374,11 +377,17 @@ def _fit_drift_params(sysid: dict, guide_exp: float, verbose: bool):
"""
free_drift_sessions = [s for s in sysid["sessions"] if s["type"] == "free_drift"]
+ # A rate from a very short window is noise, not a measurement: at a ~0.5" noise
+ # floor a 25s sample has ~±0.07 "/s slope uncertainty — one such point silently
+ # set k_ref_dec on real data (see DEC_OFFSET_ROOT_CAUSE.md).
+ MIN_FIT_SPAN_S = 120.0
+
ra_rates = []
dec_rates = []
cos2_alts = [] # 1/cos²(alt) for each session
q_factors = [] # sin(q)/cos²(alt) for each session
q_angles = []
+ fit_alts = [] # altitudes of the sessions actually used (runtime clamps to these)
for s in free_drift_sessions:
frames = s["frames"]
@@ -394,6 +403,11 @@ def _fit_drift_params(sysid: dict, guide_exp: float, verbose: bool):
print(f" [drift] skipping {s.get('session_id', '?')}: only {span:.0f}s of "
f"{requested:.0f}s requested — truncated, not a drift measurement")
continue
+ if span < MIN_FIT_SPAN_S:
+ if verbose:
+ print(f" [drift] skipping {s.get('session_id', '?')}: {span:.0f}s span — "
+ f"too short to constrain a rate (need >= {MIN_FIT_SPAN_S:.0f}s)")
+ continue
alt = s.get("altitude_deg", 45.0)
@@ -419,16 +433,21 @@ def _fit_drift_params(sysid: dict, guide_exp: float, verbose: bool):
q_angles.append(avg_q_deg)
q_rad = np.radians(avg_q_deg)
q_factors.append(np.sin(q_rad) / (cos_alt ** 2))
+ fit_alts.append(alt)
if not ra_rates:
- return 0.0, 0.0, 0.0, 0.0, 0.0
+ return 0.0, 0.0, 0.0, 0.0, 0.0, [], []
- # RA: rate = k_ref / cos²(alt) + drift_ra_extra
- if len(ra_rates) >= 2:
+ # RA: rate = k_ref / cos²(alt) + drift_ra_extra. Without real altitude spread the
+ # regression divides by ~zero variance and returns garbage — fall back to the mean.
+ alt_range = max(fit_alts) - min(fit_alts) if len(fit_alts) >= 2 else 0.0
+ if len(ra_rates) >= 2 and alt_range >= 10.0:
k_ref, drift_ra = scipy.stats.linregress(cos2_alts, ra_rates)[:2]
else:
+ if verbose and len(ra_rates) >= 2:
+ print(f" [RA] altitude range: {alt_range:.1f}° < 10° — k_ref = 0, using mean rate")
k_ref = 0.0
- drift_ra = ra_rates[0]
+ drift_ra = float(np.mean(ra_rates))
# DEC: rate = d_polar + k_ref_dec * sin(q)/cos²(alt)
q_range = max(q_angles) - min(q_angles) if len(q_angles) >= 2 else 0.0
@@ -438,18 +457,20 @@ def _fit_drift_params(sysid: dict, guide_exp: float, verbose: bool):
if verbose:
print(f" [DEC] Parallactic angle range: {q_range:.1f}° (sufficient)")
elif len(dec_rates) >= 2:
+ # No geometric spread means sin(q)/cos²(alt) is unconstrained; borrowing k_ref
+ # would apply an RA coefficient to DEC geometry it was never fitted at.
if verbose:
- print(f" [DEC] Parallactic angle range: {q_range:.1f}° < 20° — "
- f"falling back k_ref_dec = k_ref")
- k_ref_dec = k_ref
- d_polar = float(np.mean([r - k_ref * qf for r, qf in zip(dec_rates, q_factors)]))
+ print(f" [DEC] Parallactic angle range: {q_range:.1f}° < 20° — k_ref_dec = 0")
+ k_ref_dec = 0.0
+ d_polar = float(np.mean(dec_rates))
else:
d_polar = dec_rates[0] if dec_rates else 0.0
k_ref_dec = 0.0
drift_dec = 0.0 # Absorbed into d_polar
- return float(drift_ra), float(drift_dec), float(d_polar), float(k_ref), float(k_ref_dec)
+ return (float(drift_ra), float(drift_dec), float(d_polar), float(k_ref), float(k_ref_dec),
+ fit_alts, q_angles)
# ═══════════════════════════════════════════════════════════════════════════════