[education/kstars/stable-3.8.4] kstars/ekos/guide/internalguide: AI Guide: generalize harmonic PE model to N lines
Jasem Mutlaq <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit a4b5568f777342f9d7e4190bb06ffb771136978d by Jasem Mutlaq, on behalf of Pavan Kumar S G.
Committed on 15/08/2026 at 11:41.
Pushed by mutlaqja into branch 'stable-3.8.4'.
AI Guide: generalize harmonic PE model to N lines
M +110 -128 kstars/ekos/guide/internalguide/harmonic_guider.cpp
M +46 -29 kstars/ekos/guide/internalguide/harmonic_guider.h
https://invent.kde.org/education/kstars/-/commit/a4b5568f777342f9d7e4190bb06ffb771136978d
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.cpp b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
index d80e783ac2..4a493cc3ca 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.cpp
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
@@ -30,8 +30,7 @@ double HarmonicGuider::m_uncorrPosRA { 0.0 };
double HarmonicGuider::m_uncorrPosDEC { 0.0 };
int HarmonicGuider::m_frameCount { 0 };
double HarmonicGuider::m_typicalRMS { 0.5 };
-double HarmonicGuider::s_activePePeriod { -1.0 };
-double HarmonicGuider::s_activePe2Period { -1.0 };
+std::array<double, HarmonicGuider::N_LINES_MAX> HarmonicGuider::s_activeLinePeriod { };
HarmonicGuider::HarmonicGuider()
{
@@ -93,25 +92,58 @@ bool HarmonicGuider::loadWeights(const QString &weightsPath)
m_k_ref = phys["k_ref"].toDouble(0.0);
m_d_polar = phys["d_polar"].toDouble(0.0);
m_k_ref_dec = phys["k_ref_dec"].toDouble(0.0);
- m_pe_period = phys["pe_period"].toDouble(0.0);
- m_pe_amplitude = phys["pe_amplitude"].toDouble(0.0);
- // Second oscillator: the strongest secondary line the trainer found. Only accept a
- // period well separated from the primary so the two rotators cannot fight over one line.
- m_pe2_period = 0.0;
- m_pe2_amplitude = 0.0;
+ // Falls back to the legacy pe_period/pe_amplitude pair when pe_lines is absent.
+ m_linePeriod.fill(0.0);
+ m_lineAmplitude.fill(0.0);
+ m_activeLines = 0;
+
const QJsonArray peLines = phys["pe_lines"].toArray();
- if (m_pe_period > 0.0 && peLines.size() >= 2)
+ if (!peLines.isEmpty())
+ {
+ for (const QJsonValue &v : peLines)
+ {
+ if (m_activeLines >= N_LINES_MAX)
+ break;
+ const QJsonObject line = v.toObject();
+ const double period = line["period_s"].toDouble(0.0);
+ const double amplitude = line["amplitude_px"].toDouble(0.0);
+ if (period <= 0.0)
+ continue;
+ // Later lines must clear an amplitude floor and stay separated from every line already kept.
+ if (m_activeLines > 0)
+ {
+ if (amplitude <= 0.05)
+ continue;
+ bool tooClose = false;
+ for (int i = 0; i < m_activeLines; ++i)
+ {
+ if (std::abs(period - m_linePeriod[i]) / m_linePeriod[i] <= 0.15)
+ {
+ tooClose = true;
+ break;
+ }
+ }
+ if (tooClose)
+ continue;
+ }
+ m_linePeriod[m_activeLines] = period;
+ m_lineAmplitude[m_activeLines] = amplitude;
+ m_activeLines++;
+ }
+ }
+ else
{
- const QJsonObject line2 = peLines[1].toObject();
- const double p2 = line2["period_s"].toDouble(0.0);
- const double a2 = line2["amplitude_px"].toDouble(0.0);
- if (p2 > 0.0 && a2 > 0.05 && std::abs(p2 - m_pe_period) / m_pe_period > 0.15)
+ const double p0 = phys["pe_period"].toDouble(0.0);
+ const double a0 = phys["pe_amplitude"].toDouble(0.0);
+ if (p0 > 0.0)
{
- m_pe2_period = p2;
- m_pe2_amplitude = a2;
+ m_linePeriod[0] = p0;
+ m_lineAmplitude[0] = a0;
+ m_activeLines = 1;
}
}
+
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);
@@ -181,12 +213,11 @@ void HarmonicGuider::resetSession(bool forceReset)
m_qFeatRA = 0.0;
m_qFeatDec = 0.0;
- // If the loaded weights describe a different mount (different PE period), the
+ // If the loaded weights describe a different mount (different PE lines), the
// persisted static Kalman state belongs to the previous mount — discard it entirely.
- if (m_pe_period != s_activePePeriod || m_pe2_period != s_activePe2Period)
+ if (m_linePeriod != s_activeLinePeriod)
{
- s_activePePeriod = m_pe_period;
- s_activePe2Period = m_pe2_period;
+ s_activeLinePeriod = m_linePeriod;
forceReset = true;
}
@@ -237,42 +268,24 @@ void HarmonicGuider::buildF(Eigen::Matrix<double, N_STATES, N_STATES> &F, double
F(RA_SPRING, RA_SPRING) = std::exp(-dt / m_tau_ra);
F(DEC_SPRING, DEC_SPRING) = std::exp(-dt / m_tau_dec);
- // PE states evolve as 2D rotation at the PE frequency
- if (m_pe_period > 0.0)
+ // Each active PE line evolves as its own 2D rotation, independent of the others.
+ for (int l = 0; l < m_activeLines; ++l)
{
- const double omega = 2.0 * M_PI / m_pe_period;
- const double cos_wdt = std::cos(omega * dt);
- const double sin_wdt = std::sin(omega * dt);
-
- // RA PE rotation
- F(RA_PE_SIN, RA_PE_SIN) = cos_wdt;
- F(RA_PE_SIN, RA_PE_COS) = sin_wdt;
- F(RA_PE_COS, RA_PE_SIN) = -sin_wdt;
- F(RA_PE_COS, RA_PE_COS) = cos_wdt;
-
- // DEC PE rotation
- F(DEC_PE_SIN, DEC_PE_SIN) = cos_wdt;
- F(DEC_PE_SIN, DEC_PE_COS) = sin_wdt;
- F(DEC_PE_COS, DEC_PE_SIN) = -sin_wdt;
- F(DEC_PE_COS, DEC_PE_COS) = cos_wdt;
- }
+ const double omega = 2.0 * M_PI / m_linePeriod[l];
+ const double c = std::cos(omega * dt);
+ const double s = std::sin(omega * dt);
+ const int raS = raSinIdx(l), raC = raCosIdx(l);
+ const int deS = decSinIdx(l), deC = decCosIdx(l);
- // Second PE line rotates at its own frequency
- if (m_pe2_period > 0.0)
- {
- const double omega2 = 2.0 * M_PI / m_pe2_period;
- const double c2 = std::cos(omega2 * dt);
- const double s2 = std::sin(omega2 * dt);
-
- F(RA_PE2_SIN, RA_PE2_SIN) = c2;
- F(RA_PE2_SIN, RA_PE2_COS) = s2;
- F(RA_PE2_COS, RA_PE2_SIN) = -s2;
- F(RA_PE2_COS, RA_PE2_COS) = c2;
-
- F(DEC_PE2_SIN, DEC_PE2_SIN) = c2;
- F(DEC_PE2_SIN, DEC_PE2_COS) = s2;
- F(DEC_PE2_COS, DEC_PE2_SIN) = -s2;
- F(DEC_PE2_COS, DEC_PE2_COS) = c2;
+ F(raS, raS) = c;
+ F(raS, raC) = s;
+ F(raC, raS) = -s;
+ F(raC, raC) = c;
+
+ F(deS, deS) = c;
+ F(deS, deC) = s;
+ F(deC, deS) = -s;
+ F(deC, deC) = c;
}
}
@@ -315,20 +328,14 @@ 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)
- if (m_pe_period > 0.0)
+ // PE process noise (small — PE is nearly deterministic), same fixed constant
+ // for every active line; per-line scaling is a possible future refinement.
+ for (int l = 0; l < m_activeLines; ++l)
{
- 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;
+ Q(raSinIdx(l), raSinIdx(l)) = 0.001 * dt;
+ Q(raCosIdx(l), raCosIdx(l)) = 0.001 * dt;
+ Q(decSinIdx(l), decSinIdx(l)) = 0.001 * dt;
+ Q(decCosIdx(l), decCosIdx(l)) = 0.001 * dt;
}
return Q;
@@ -354,6 +361,22 @@ void HarmonicGuider::driftRates(double alt_deg, double parallactic_angle_deg,
}
}
+double HarmonicGuider::peSumRA(const Eigen::Matrix<double, N_STATES, 1> &x) const
+{
+ double sum = 0.0;
+ for (int l = 0; l < m_activeLines; ++l)
+ sum += x(raSinIdx(l));
+ return sum;
+}
+
+double HarmonicGuider::peSumDEC(const Eigen::Matrix<double, N_STATES, 1> &x) const
+{
+ double sum = 0.0;
+ for (int l = 0; l < m_activeLines; ++l)
+ sum += x(decSinIdx(l));
+ return sum;
+}
+
// ── Kalman predict step: free dynamics only, pulses are applied in update() ──
void HarmonicGuider::kalmanPredict(double dt, double alt_deg, double parallactic_angle_deg)
{
@@ -390,20 +413,15 @@ 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
- // H extracts: ra_obs = ra_err + pe_sin_ra, dec_obs = dec_err + pe_sin_dec
+ // Observation matrix: observe position + the sin component of every active PE line.
+ // ra_obs = ra_err + Σ pe_sin_ra[l], dec_obs = dec_err + Σ pe_sin_dec[l]
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)
+ for (int l = 0; l < m_activeLines; ++l)
{
- H(0, RA_PE2_SIN) = 1.0;
- H(1, DEC_PE2_SIN) = 1.0;
+ H(0, raSinIdx(l)) = 1.0;
+ H(1, decSinIdx(l)) = 1.0;
}
// Measurement noise from the current frame's SNR (~0.5 px at SNR 30)
@@ -418,7 +436,7 @@ void HarmonicGuider::kalmanUpdate(double ra_meas_px, double dec_meas_px, double
// Innovation covariance
Eigen::Matrix2d S = H * m_P * H.transpose() + R;
- // Kalman gain (2x10 matrix)
+ // Kalman gain
Eigen::Matrix<double, N_STATES, N_OBS> K = m_P * H.transpose() * S.inverse();
// State update
@@ -471,33 +489,13 @@ GuideOutput HarmonicGuider::predict(const GuideFrameData &frame)
m_hasPrevRaw = true;
// Posterior (POS + PE) before propagation
- double post_ra = m_x(RA_POS);
- double post_dec = m_x(DEC_POS);
- if (m_pe_period > 0.0)
- {
- post_ra += m_x(RA_PE_SIN);
- post_dec += m_x(DEC_PE_SIN);
- }
- if (m_pe2_period > 0.0)
- {
- post_ra += m_x(RA_PE2_SIN);
- post_dec += m_x(DEC_PE2_SIN);
- }
+ const double post_ra = m_x(RA_POS) + peSumRA(m_x);
+ const double post_dec = m_x(DEC_POS) + peSumDEC(m_x);
kalmanPredict(frame.dt, frame.altitude_deg, frame.parallactic_angle_deg);
- double pred_ra = m_x(RA_POS);
- double pred_dec = m_x(DEC_POS);
- if (m_pe_period > 0.0)
- {
- pred_ra += m_x(RA_PE_SIN);
- pred_dec += m_x(DEC_PE_SIN);
- }
- if (m_pe2_period > 0.0)
- {
- pred_ra += m_x(RA_PE2_SIN);
- pred_dec += m_x(DEC_PE2_SIN);
- }
+ const double pred_ra = m_x(RA_POS) + peSumRA(m_x);
+ const double pred_dec = m_x(DEC_POS) + peSumDEC(m_x);
// Prediction is the expected uncorrected drift over the next interval
m_lastPredRA = pred_ra - post_ra;
@@ -558,9 +556,13 @@ void HarmonicGuider::update(double /*ra_error_px*/, double /*dec_error_px*/,
QString HarmonicGuider::stateString() const
{
- QString pe = QString::number(m_pe_period, 'f', 1);
- if (m_pe2_period > 0.0)
- pe += QString("+%1").arg(m_pe2_period, 0, 'f', 1);
+ QString pe;
+ for (int l = 0; l < m_activeLines; ++l)
+ {
+ if (l > 0) pe += "+";
+ pe += QString::number(m_linePeriod[l], 'f', 1);
+ }
+ if (pe.isEmpty()) pe = "0";
return QString("Harmonic κ_ra=%1 τ_ra=%2 PE=%3s conf=%4")
.arg(m_kappa_ra, 0, 'f', 2)
.arg(m_tau_ra, 0, 'f', 1)
@@ -623,18 +625,8 @@ GuideOutput HarmonicGuider::darkPredict(double dt_sec)
Eigen::Matrix<double, N_STATES, N_STATES> F;
buildF(F, dt_sec);
- double post_ra = x(RA_POS);
- double post_dec = x(DEC_POS);
- if (m_pe_period > 0.0)
- {
- post_ra += x(RA_PE_SIN);
- post_dec += x(DEC_PE_SIN);
- }
- if (m_pe2_period > 0.0)
- {
- post_ra += x(RA_PE2_SIN);
- post_dec += x(DEC_PE2_SIN);
- }
+ const double post_ra = x(RA_POS) + peSumRA(x);
+ const double post_dec = x(DEC_POS) + peSumDEC(x);
const double release_ra = x(RA_SPRING) * (1.0 - std::exp(-dt_sec / m_tau_ra));
const double release_dec = x(DEC_SPRING) * (1.0 - std::exp(-dt_sec / m_tau_dec));
@@ -647,18 +639,8 @@ GuideOutput HarmonicGuider::darkPredict(double dt_sec)
x(RA_POS) += ra_rate * dt_sec;
x(DEC_POS) += dec_rate * dt_sec;
- double pred_ra = x(RA_POS);
- double pred_dec = x(DEC_POS);
- if (m_pe_period > 0.0)
- {
- pred_ra += x(RA_PE_SIN);
- pred_dec += x(DEC_PE_SIN);
- }
- if (m_pe2_period > 0.0)
- {
- pred_ra += x(RA_PE2_SIN);
- pred_dec += x(DEC_PE2_SIN);
- }
+ const double pred_ra = x(RA_POS) + peSumRA(x);
+ const double pred_dec = x(DEC_POS) + peSumDEC(x);
GuideOutput out;
out.valid = (m_frameCount > warmupFrames());
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.h b/kstars/ekos/guide/internalguide/harmonic_guider.h
index a7b3f71583..01ba47591a 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.h
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.h
@@ -2,24 +2,26 @@
/*
* harmonic_guider.h — Neural Kalman Filter guider for harmonic drive mounts
*
- * Architecture: 14-state Kalman filter with spring dynamics and two PE oscillators,
- * plus a small Q-net MLP (~66 parameters) for adaptive process noise.
+ * Architecture: Kalman filter with spring dynamics and up to N_LINES_MAX PE
+ * oscillators per axis, plus a small Q-net MLP (~66 parameters)
+ * for adaptive process noise.
*
- * State vector:
- * [ra_err, ra_vel, spring_ra, pe_sin_ra, pe_cos_ra,
- * dec_err, dec_vel, spring_dec, pe_sin_dec, pe_cos_dec,
- * pe2_sin_ra, pe2_cos_ra, pe2_sin_dec, pe2_cos_dec]
+ * State vector: [ra_err, ra_vel, ra_spring, dec_err, dec_vel, dec_spring,
+ * (pe_sin, pe_cos) x N_LINES_MAX for RA,
+ * (pe_sin, pe_cos) x N_LINES_MAX for DEC]
+ * An inactive line (period == 0) contributes identity/zero to F/Q/H.
*
* The spring states model elastic wind-up: the flexspline absorbs a fraction κ
* of each correction pulse and releases it exponentially with time constant τ.
*
- * 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.
+ * Each PE line is an independent 2D rotation at its own frequency; lines need
+ * not be harmonically related.
*
* SPDX-License-Identifier: GPL-2.0-or-later
*/
#include "mount_guider.h"
+#include <array>
#include <deque>
#include <Eigen/Core>
#include <QJsonObject>
@@ -53,24 +55,37 @@ class HarmonicGuider : public MountSpecificGuider
private:
// ── Kalman dimensions ────────────────────────────────────────────────
- static constexpr int N_STATES = 14;
+ /// Cap on simultaneously-tracked PE lines; matches the trainer's cap
+ /// (train_harmonic.py::_estimate_pe keeps at most 4 surviving lines).
+ static constexpr int N_LINES_MAX = 4;
+ static constexpr int N_STATES = 6 + 4 * N_LINES_MAX;
static constexpr int N_OBS = 2;
- // State indices
+ // Fixed state indices
static constexpr int RA_POS = 0;
static constexpr int RA_VEL = 1;
static constexpr int RA_SPRING = 2;
- static constexpr int RA_PE_SIN = 3;
- static constexpr int RA_PE_COS = 4;
- static constexpr int DEC_POS = 5;
- static constexpr int DEC_VEL = 6;
- static constexpr int DEC_SPRING = 7;
- static constexpr int DEC_PE_SIN = 8;
- static constexpr int DEC_PE_COS = 9;
- static constexpr int RA_PE2_SIN = 10;
- static constexpr int RA_PE2_COS = 11;
- static constexpr int DEC_PE2_SIN = 12;
- static constexpr int DEC_PE2_COS = 13;
+ static constexpr int DEC_POS = 3;
+ static constexpr int DEC_VEL = 4;
+ static constexpr int DEC_SPRING = 5;
+
+ // PE line l (0-based) occupies four states from this base offset.
+ static constexpr int raSinIdx(int l)
+ {
+ return 6 + 4 * l;
+ }
+ static constexpr int raCosIdx(int l)
+ {
+ return 7 + 4 * l;
+ }
+ static constexpr int decSinIdx(int l)
+ {
+ return 8 + 4 * l;
+ }
+ static constexpr int decCosIdx(int l)
+ {
+ return 9 + 4 * l;
+ }
// ── Spring parameters (loaded from weights JSON) ─────────────────────
double m_kappa_ra { 0.2 }; ///< Fraction of RA pulse absorbed by spring [0, 0.9]
@@ -78,11 +93,11 @@ class HarmonicGuider : public MountSpecificGuider
double m_kappa_dec { 0.2 }; ///< Fraction of DEC pulse absorbed
double m_tau_dec { 1.5 }; ///< DEC spring release time constant
- // ── PE parameters (loaded from weights JSON) ─────────────────────────
- double m_pe_period { 0.0 }; ///< PE period in seconds (0 = no PE detected)
- double m_pe_amplitude { 0.0 }; ///< Initial PE amplitude estimate (px)
- double m_pe2_period { 0.0 }; ///< Second PE line period (0 = single-line model)
- double m_pe2_amplitude { 0.0 };
+ // ── PE parameters (loaded from weights JSON) ──────────────────────────
+ // Slot 0 is the primary line; a period of 0 means that slot is inactive.
+ std::array<double, N_LINES_MAX> m_linePeriod { };
+ std::array<double, N_LINES_MAX> m_lineAmplitude { };
+ int m_activeLines { 0 }; ///< count of slots with period > 0, always a prefix
// ── Drift / refraction parameters ────────────────────────────────────
double m_drift_ra { 0.0 }; ///< Baseline RA drift rate (px/s)
@@ -102,10 +117,9 @@ class HarmonicGuider : public MountSpecificGuider
static double m_uncorrPosDEC;
static int m_frameCount;
static double m_typicalRMS;
- /// PE period the persisted static state was built for; a change means the
+ /// Line periods the persisted static state was built for; a change means the
/// loaded weights describe a different mount, so the static state is discarded.
- static double s_activePePeriod;
- static double s_activePe2Period;
+ static std::array<double, N_LINES_MAX> s_activeLinePeriod;
// ── Q-net MLP weights (5 → 8 → 2) ───────────────────────────────────
// Input: [snr_norm, snr_delta_norm, |innov_ra|, |innov_dec|, dt_norm]
@@ -156,4 +170,7 @@ class HarmonicGuider : public MountSpecificGuider
void kalmanPredict(double dt, double alt_deg, double parallactic_angle_deg);
void kalmanUpdate(double ra_meas_px, double dec_meas_px, double snr);
void updateConfidence(double innovRA, double innovDec, double snr);
+ /// Sum of the sin-component (position contribution) of every active PE line.
+ double peSumRA(const Eigen::Matrix<double, N_STATES, 1> &x) const;
+ double peSumDEC(const Eigen::Matrix<double, N_STATES, 1> &x) const;
};