[education/kstars] kstars/ekos/guide: AI Guide: Fix filter models, rework sysid protocol and trainer, improve fingerprint diagnostics
Jasem Mutlaq <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit e65e53472adb3c98a053f58fb053f162fd450568 by Jasem Mutlaq, on behalf of Pavan Kumar S G.
Committed on 23/07/2026 at 18:32.
Pushed by mutlaqja into branch 'master'.
AI Guide: Fix filter models, rework sysid protocol and trainer, improve fingerprint diagnostics
M +49 -19 kstars/ekos/guide/aiguideprotocol.cpp
M +2 -0 kstars/ekos/guide/aiguideprotocol.h
M +1 -1 kstars/ekos/guide/aiguidewizard.ui
M +34 -19 kstars/ekos/guide/internalguide/direct_drive_guider.cpp
M +4 -2 kstars/ekos/guide/internalguide/direct_drive_guider.h
M +3 -1 kstars/ekos/guide/internalguide/gmath.cpp
M +139 -142 kstars/ekos/guide/internalguide/harmonic_guider.cpp
M +10 -4 kstars/ekos/guide/internalguide/harmonic_guider.h
M +4 -1 kstars/ekos/guide/internalguide/mount_guider.h
M +42 -50 kstars/ekos/guide/internalguide/worm_gear_guider.cpp
M +3 -1 kstars/ekos/guide/internalguide/worm_gear_guider.h
M +2 -0 kstars/ekos/guide/offlinetrainer/train_direct_drive.py
M +259 -133 kstars/ekos/guide/offlinetrainer/train_harmonic.py
M +5 -1 kstars/ekos/guide/offlinetrainer/train_worm_gear.py
https://invent.kde.org/education/kstars/-/commit/e65e53472adb3c98a053f58fb053f162fd450568
diff --git a/kstars/ekos/guide/aiguideprotocol.cpp b/kstars/ekos/guide/aiguideprotocol.cpp
index a692c5b9c6..5c42ec7ca6 100644
--- a/kstars/ekos/guide/aiguideprotocol.cpp
+++ b/kstars/ekos/guide/aiguideprotocol.cpp
@@ -126,8 +126,8 @@ void AIGuideProtocol::start(const QString &mountType)
fingerprint["dec_integral_gain"] = Options::dECIntegralGain();
fingerprint["ra_min_pulse_arcsec"] = Options::rAMinimumPulseArcSec();
fingerprint["dec_min_pulse_arcsec"] = Options::dECMinimumPulseArcSec();
- fingerprint["ra_max_pulse_arcsec"] = static_cast<int>(Options::rAMaximumPulseArcSec());
- fingerprint["dec_max_pulse_arcsec"] = static_cast<int>(Options::dECMaximumPulseArcSec());
+ fingerprint["ra_max_pulse_arcsec"] = static_cast<double>(Options::rAMaximumPulseArcSec());
+ fingerprint["dec_max_pulse_arcsec"] = static_cast<double>(Options::dECMaximumPulseArcSec());
fingerprint["ra_hysteresis"] = Options::rAHysteresis();
fingerprint["dec_hysteresis"] = Options::dECHysteresis();
fingerprint["ra_pulse_algorithm"] = 0;
@@ -154,22 +154,26 @@ void AIGuideProtocol::start(const QString &mountType)
}
else if (mountStr == "Harmonic Drive")
{
- m_Phases.append({65.0, -45.0, 480, false, false, {}, {}, 0, 15, 30});
- m_Phases.append({65.0, -45.0, 120, true, false, {}, {}, 0, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "EAST", 50, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "EAST", 100, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "EAST", 200, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "WEST", 50, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "WEST", 100, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "RA", "WEST", 200, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "NORTH", 50, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "NORTH", 100, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "NORTH", 200, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "SOUTH", 50, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "SOUTH", 100, 15, 30});
- m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "SOUTH", 200, 15, 30});
- m_Phases.append({45.0, -45.0, 120, true, false, {}, {}, 0, 15, 30});
- m_Phases.append({45.0, -45.0, 300, false, false, {}, {}, 0, 15, 30});
+ // Position 1: 720s standard guiding and 480 free drift
+ m_Phases.append({65.0, -45.0, 720, false, false, {}, {}, 0, 15, 20});
+ m_Phases.append({65.0, -45.0, 480, true, false, {}, {}, 0, 15, 20});
+
+ // Pulse response: large pulses, alternating direction so drift/PE cancel in pairing
+ for (int rep = 0; rep < 3; rep++)
+ {
+ m_Phases.append({65.0, -45.0, 0, false, true, "RA", "EAST", 500, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "RA", "WEST", 500, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "RA", "EAST", 1000, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "RA", "WEST", 1000, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "NORTH", 500, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "SOUTH", 500, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "NORTH", 1000, 15, 20});
+ m_Phases.append({65.0, -45.0, 0, false, true, "DEC", "SOUTH", 1000, 15, 20});
+ }
+
+ // Position 2 east of the meridian: parallactic spread for the DEC refraction fit
+ m_Phases.append({45.0, 45.0, 120, true, false, {}, {}, 0, 15, 20});
+ m_Phases.append({45.0, 45.0, 300, false, false, {}, {}, 0, 15, 20});
}
else
{
@@ -569,6 +573,7 @@ void AIGuideProtocol::processProtocol()
m_Guide->setAIFreeDrift(true);
m_PulseFrameCount = 0;
m_PulseResponseData = QJsonArray();
+ m_PulseBaselineData = QJsonArray();
connect(m_Guide, &Guide::guideStats, this, &AIGuideProtocol::onGuideStats, Qt::UniqueConnection);
@@ -614,6 +619,7 @@ void AIGuideProtocol::processProtocol()
m_PulseFrameCount = 0;
m_PulseResponseData = QJsonArray();
+ m_PulseSentAtMs = QDateTime::currentMSecsSinceEpoch();
m_FrameTimer.start();
m_PulseWatchdog = phase.responseFrames * 6 + 30;
m_State = STATE_PULSE_RECORDING;
@@ -663,6 +669,8 @@ void AIGuideProtocol::processProtocol()
pulseSession["pier_side"] = (pierSide == ISD::Mount::PIER_EAST) ? "EAST" : "WEST";
}
+ pulseSession["baseline_frames"] = m_PulseBaselineData;
+ pulseSession["pulse_sent_at_ms"] = m_PulseSentAtMs;
pulseSession["response_frames"] = m_PulseResponseData;
QJsonArray sessions = m_SysIdData["sessions"].toArray();
@@ -802,6 +810,27 @@ void AIGuideProtocol::onGuideStats(double raErr, double decErr, int raPulse, int
m_PhaseData.append(frame);
}
+ // Pre-pulse settle frames become the fit's baseline
+ if (m_State == STATE_PULSE_SETTLING && m_PulseFrameCount == 0 && m_PulseResponseData.isEmpty()
+ && m_PulseBaselineData.size() < 6)
+ {
+ double dx = raErr;
+ double dy = decErr;
+ if (m_Guide && m_Guide->pixelSizeX() > 0 && m_Guide->focalLength() > 0)
+ {
+ double binning = std::max(1, Options::guideBinning().left(1).toInt());
+ double scale = (206.265 * m_Guide->pixelSizeX() * binning) / m_Guide->focalLength();
+ dx = -raErr / scale;
+ dy = decErr / scale;
+ }
+
+ QJsonObject frame;
+ frame["ra_raw_px"] = dx;
+ frame["dec_raw_px"] = dy;
+ frame["snr"] = snr;
+ m_PulseBaselineData.append(frame);
+ }
+
if (m_State == STATE_PULSE_RECORDING)
{
double dt = m_FrameTimer.isValid() ? (m_FrameTimer.restart() / 1000.0) : 0.0;
@@ -817,7 +846,8 @@ void AIGuideProtocol::onGuideStats(double raErr, double decErr, int raPulse, int
}
QJsonObject frame;
- frame["t"] = m_PulseFrameCount * dt;
+ // True seconds since the pulse was sent
+ frame["t"] = (QDateTime::currentMSecsSinceEpoch() - m_PulseSentAtMs) / 1000.0;
frame["ra_raw_px"] = dx;
frame["dec_raw_px"] = dy;
frame["snr"] = snr;
diff --git a/kstars/ekos/guide/aiguideprotocol.h b/kstars/ekos/guide/aiguideprotocol.h
index 51cdf68188..e462bff13d 100644
--- a/kstars/ekos/guide/aiguideprotocol.h
+++ b/kstars/ekos/guide/aiguideprotocol.h
@@ -151,6 +151,8 @@ class AIGuideProtocol : public QObject
int m_PulseSettleTimer { 0 };
int m_PulseWatchdog { 0 };
QJsonArray m_PulseResponseData;
+ QJsonArray m_PulseBaselineData; ///< pre-pulse frames, the fit's reference level
+ qint64 m_PulseSentAtMs { 0 }; ///< pulse send time, t=0 for response frames
};
}
\ No newline at end of file
diff --git a/kstars/ekos/guide/aiguidewizard.ui b/kstars/ekos/guide/aiguidewizard.ui
index 51cee1bea3..ccf03bff98 100644
--- a/kstars/ekos/guide/aiguidewizard.ui
+++ b/kstars/ekos/guide/aiguidewizard.ui
@@ -74,7 +74,7 @@
<item>
<widget class="QLabel" name="exposureLabel">
<property name="text">
[suppressed due to size limit]
[suppressed due to size limit]
</property>
<property name="wordWrap">
<bool>true</bool>
diff --git a/kstars/ekos/guide/internalguide/direct_drive_guider.cpp b/kstars/ekos/guide/internalguide/direct_drive_guider.cpp
index 734fd5380c..ada4f3fdee 100644
--- a/kstars/ekos/guide/internalguide/direct_drive_guider.cpp
+++ b/kstars/ekos/guide/internalguide/direct_drive_guider.cpp
@@ -7,6 +7,7 @@
#include "direct_drive_guider.h"
#include "Options.h"
+#include "ekos_guide_debug.h"
#include <QFile>
#include <QJsonDocument>
@@ -31,26 +32,37 @@ bool DirectDriveGuider::validateFingerprint(const QJsonObject &fp)
if (fp.isEmpty())
return true;
- if (fp.contains("guide_exposure_s"))
+ const struct
{
- const double expected = fp["guide_exposure_s"].toDouble();
- if (!fpDoubleClose(expected, Options::guideExposure(), 0.05))
- return false;
- }
-
- if (fp.contains("guide_binning") &&
- fp["guide_binning"].toString() != Options::guideBinning())
- return false;
+ const char *key;
+ double current;
+ double tol;
+ } checks[] =
+ {
+ { "guide_exposure_s", Options::guideExposure(), 0.05 },
+ { "ra_proportional_gain", Options::rAProportionalGain(), 1e-4 },
+ { "dec_proportional_gain", Options::dECProportionalGain(), 1e-4 },
+ };
- if (fp.contains("ra_proportional_gain") &&
- !fpDoubleClose(fp["ra_proportional_gain"].toDouble(), Options::rAProportionalGain()))
- return false;
+ bool ok = true;
+ for (const auto &c : checks)
+ {
+ if (fp.contains(c.key) && !fpDoubleClose(fp[c.key].toDouble(), c.current, c.tol))
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected:" << c.key << "recorded"
+ << fp[c.key].toDouble() << "current" << c.current;
+ ok = false;
+ }
+ }
- if (fp.contains("dec_proportional_gain") &&
- !fpDoubleClose(fp["dec_proportional_gain"].toDouble(), Options::dECProportionalGain()))
- return false;
+ if (fp.contains("guide_binning") && fp["guide_binning"].toString() != Options::guideBinning())
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected: guide_binning recorded"
+ << fp["guide_binning"].toString() << "current" << Options::guideBinning();
+ ok = false;
+ }
- return true;
+ return ok;
}
// ---------------------------------------------------------------------------
@@ -85,6 +97,8 @@ bool DirectDriveGuider::loadWeights(const QString &weightsPath)
m_k_ref_dec = params["k_ref_dec"].toDouble(0.0);
m_d_ra_extra = params["d_ra_extra"].toDouble(0.0);
m_phi_drift = params["phi_drift"].toDouble(0.0);
+ m_fit_alt_min = params["fit_alt_min"].toDouble(35.0);
+ m_fit_alt_max = params["fit_alt_max"].toDouble(65.0);
m_pixel_scale = root["pixel_scale"].toDouble(1.0);
@@ -95,7 +109,7 @@ bool DirectDriveGuider::loadWeights(const QString &weightsPath)
// ---------------------------------------------------------------------------
// resetSession
// ---------------------------------------------------------------------------
-void DirectDriveGuider::resetSession(bool forceReset)
+void DirectDriveGuider::resetSession(bool /*forceReset*/)
{
m_frameCount = 0;
}
@@ -169,7 +183,8 @@ GuideOutput DirectDriveGuider::darkPredict(double dt_sec)
// ---------------------------------------------------------------------------
double DirectDriveGuider::physicsRA(double alt_deg) const
{
- const double alt_rad = alt_deg * M_PI / 180.0;
+ // Refraction fit is only valid inside the fitted altitude range.
+ 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);
if (std::abs(cos_alt) < 1e-4) return m_d_ra_extra;
return m_k_ref / (cos_alt * cos_alt) + m_d_ra_extra;
@@ -180,7 +195,7 @@ double DirectDriveGuider::physicsRA(double alt_deg) const
// ---------------------------------------------------------------------------
double DirectDriveGuider::physicsDEC(double alt_deg, double q_deg) const
{
- const double alt_rad = alt_deg * M_PI / 180.0;
+ const double alt_rad = std::clamp(alt_deg, m_fit_alt_min, m_fit_alt_max) * M_PI / 180.0;
const double q_rad = q_deg * M_PI / 180.0;
const double cos_alt = std::cos(alt_rad);
if (std::abs(cos_alt) < 1e-4) return m_d_polar;
diff --git a/kstars/ekos/guide/internalguide/direct_drive_guider.h b/kstars/ekos/guide/internalguide/direct_drive_guider.h
index 0fb4d21b46..6f6936ecca 100644
--- a/kstars/ekos/guide/internalguide/direct_drive_guider.h
+++ b/kstars/ekos/guide/internalguide/direct_drive_guider.h
@@ -28,7 +28,8 @@ class DirectDriveGuider : public MountSpecificGuider
// DirectDriveGuider update() is a no-op: there are no learnable online parameters.
void update(double /*ra_px*/, double /*dec_px*/,
- double /*uncorr_ra*/, double /*uncorr_dec*/, double /*snr*/) override {}
+ double /*uncorr_ra*/, double /*uncorr_dec*/, double /*snr*/,
+ double /*ra_pulse_px*/, double /*dec_pulse_px*/) override {}
double confidence() const override
{
@@ -42,7 +43,6 @@ class DirectDriveGuider : public MountSpecificGuider
{
return m_weightsLoaded;
}
-
QString stateString() const override;
private:
@@ -57,6 +57,8 @@ class DirectDriveGuider : public MountSpecificGuider
double m_d_ra_extra { 0.0 }; ///< Residual RA drift not explained by refraction (px/s)
double m_d_polar { 0.0 }; ///< Constant DEC polar drift (px/s)
double m_k_ref_dec { 0.0 }; ///< DEC refraction coefficient
+ double m_fit_alt_min { 35.0 }; ///< Altitude range the fit is valid for
+ double m_fit_alt_max { 65.0 };
double m_phi_drift { 0.0 }; ///< Polar drift vector angle (unused in inference, stored for info)
double m_pixel_scale { 1.0 }; ///< arcsec/pixel, loaded from weights JSON
diff --git a/kstars/ekos/guide/internalguide/gmath.cpp b/kstars/ekos/guide/internalguide/gmath.cpp
index 9b5a047fa6..f1bad9e18b 100644
--- a/kstars/ekos/guide/internalguide/gmath.cpp
+++ b/kstars/ekos/guide/internalguide/gmath.cpp
@@ -1049,7 +1049,9 @@ void cgmath::performProcessing(Ekos::GuideState state, QSharedPointer<FITSData>
double uncorrected_drift_ra_px = uncorrected_drift_ra_arcsec / frameData.pixel_scale;
double uncorrected_drift_dec_px = uncorrected_drift_dec_arcsec / frameData.pixel_scale;
- m_AIGuider->update(ra_px, dec_px, uncorrected_drift_ra_px, uncorrected_drift_dec_px, frameData.snr);
+ m_AIGuider->update(ra_px, dec_px, uncorrected_drift_ra_px, uncorrected_drift_dec_px, frameData.snr,
+ applied_pulse_arcsec_ra / frameData.pixel_scale,
+ applied_pulse_arcsec_dec / frameData.pixel_scale);
m_lastAIPrediction = m_AIGuider->predict(frameData);
qCDebug(KSTARS_EKOS_GUIDE) <<
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.cpp b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
index 1b45fb77cc..70fd9f4423 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.cpp
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.cpp
@@ -7,6 +7,7 @@
#include "harmonic_guider.h"
#include "Options.h"
+#include "ekos_guide_debug.h"
#include <Eigen/LU>
@@ -33,6 +34,8 @@ HarmonicGuider::m_x { Eigen::Matrix<double, N_STATES, 1>::Zero() };
Eigen::Matrix<double, HarmonicGuider::N_STATES, HarmonicGuider::N_STATES>
HarmonicGuider::m_P { Eigen::Matrix<double, N_STATES, N_STATES>::Identity() * 10.0 };
+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 };
@@ -55,58 +58,45 @@ bool HarmonicGuider::validateFingerprint(const QJsonObject &fp)
if (fp.isEmpty())
return true;
- if (fp.contains("guide_exposure_s"))
+ const struct
{
- const double expected = fp["guide_exposure_s"].toDouble();
- if (!fpDoubleClose(expected, Options::guideExposure(), 0.05))
- return false;
- }
-
- if (fp.contains("ra_proportional_gain") &&
- !fpDoubleClose(fp["ra_proportional_gain"].toDouble(), Options::rAProportionalGain()))
- return false;
-
- if (fp.contains("dec_proportional_gain") &&
- !fpDoubleClose(fp["dec_proportional_gain"].toDouble(), Options::dECProportionalGain()))
- return false;
-
- if (fp.contains("ra_integral_gain") &&
- !fpDoubleClose(fp["ra_integral_gain"].toDouble(), Options::rAIntegralGain()))
- return false;
-
- if (fp.contains("dec_integral_gain") &&
- !fpDoubleClose(fp["dec_integral_gain"].toDouble(), Options::dECIntegralGain()))
- return false;
-
- if (fp.contains("ra_min_pulse_arcsec") &&
- !fpDoubleClose(fp["ra_min_pulse_arcsec"].toDouble(), Options::rAMinimumPulseArcSec()))
- return false;
-
- if (fp.contains("dec_min_pulse_arcsec") &&
- !fpDoubleClose(fp["dec_min_pulse_arcsec"].toDouble(), Options::dECMinimumPulseArcSec()))
- return false;
-
- if (fp.contains("ra_max_pulse_arcsec") &&
- !fpDoubleClose(fp["ra_max_pulse_arcsec"].toDouble(), Options::rAMaximumPulseArcSec()))
- return false;
-
- if (fp.contains("dec_max_pulse_arcsec") &&
- !fpDoubleClose(fp["dec_max_pulse_arcsec"].toDouble(), Options::dECMaximumPulseArcSec()))
- return false;
-
- if (fp.contains("ra_hysteresis") &&
- !fpDoubleClose(fp["ra_hysteresis"].toDouble(), Options::rAHysteresis()))
- return false;
+ const char *key;
+ double current;
+ double tol;
+ } checks[] =
+ {
+ { "guide_exposure_s", Options::guideExposure(), 0.05 },
+ { "ra_proportional_gain", Options::rAProportionalGain(), 1e-4 },
+ { "dec_proportional_gain", Options::dECProportionalGain(), 1e-4 },
+ { "ra_integral_gain", Options::rAIntegralGain(), 1e-4 },
+ { "dec_integral_gain", Options::dECIntegralGain(), 1e-4 },
+ { "ra_min_pulse_arcsec", Options::rAMinimumPulseArcSec(), 1e-4 },
+ { "dec_min_pulse_arcsec", Options::dECMinimumPulseArcSec(), 1e-4 },
+ { "ra_max_pulse_arcsec", static_cast<double>(Options::rAMaximumPulseArcSec()), 1e-4 },
+ { "dec_max_pulse_arcsec", static_cast<double>(Options::dECMaximumPulseArcSec()), 1e-4 },
+ { "ra_hysteresis", Options::rAHysteresis(), 1e-4 },
+ { "dec_hysteresis", Options::dECHysteresis(), 1e-4 },
+ };
- if (fp.contains("dec_hysteresis") &&
- !fpDoubleClose(fp["dec_hysteresis"].toDouble(), Options::dECHysteresis()))
- return false;
+ bool ok = true;
+ for (const auto &c : checks)
+ {
+ if (fp.contains(c.key) && !fpDoubleClose(fp[c.key].toDouble(), c.current, c.tol))
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected:" << c.key << "recorded"
+ << fp[c.key].toDouble() << "current" << c.current;
+ ok = false;
+ }
+ }
- if (fp.contains("guide_binning") &&
- fp["guide_binning"].toString() != Options::guideBinning())
- return false;
+ if (fp.contains("guide_binning") && fp["guide_binning"].toString() != Options::guideBinning())
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected: guide_binning recorded"
+ << fp["guide_binning"].toString() << "current" << Options::guideBinning();
+ ok = false;
+ }
- return true;
+ return ok;
}
bool HarmonicGuider::loadWeights(const QString &weightsPath)
@@ -141,6 +131,8 @@ bool HarmonicGuider::loadWeights(const QString &weightsPath)
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);
+ m_fit_alt_min = phys["fit_alt_min"].toDouble(35.0);
+ m_fit_alt_max = phys["fit_alt_max"].toDouble(65.0);
// Sanity bounds on spring parameters
m_kappa_ra = std::clamp(m_kappa_ra, 0.0, 0.9);
@@ -214,6 +206,9 @@ void HarmonicGuider::resetSession(bool forceReset)
forceReset = true;
}
+ m_uncorrPosRA = 0.0;
+ m_uncorrPosDEC = 0.0;
+
if (forceReset)
{
// Full reset: first use, meridian flip, or mount/weights change.
@@ -310,9 +305,9 @@ HarmonicGuider::computeQ(double snr, double snr_delta,
Q(RA_POS, RA_POS) = q_ra * dt;
Q(DEC_POS, DEC_POS) = q_dec * dt;
- // Velocity process noise (smaller — velocity changes slowly)
- Q(RA_VEL, RA_VEL) = q_ra * 0.01 * dt;
- Q(DEC_VEL, DEC_VEL) = q_dec * 0.01 * dt;
+ // Fixed small noise: velocity tracks only the slow residual drift trend
+ Q(RA_VEL, RA_VEL) = 1e-5 * dt;
+ Q(DEC_VEL, DEC_VEL) = 1e-5 * dt;
// Spring process noise (very small — spring params are stable)
Q(RA_SPRING, RA_SPRING) = 0.001 * dt;
@@ -330,67 +325,43 @@ HarmonicGuider::computeQ(double snr, double snr_delta,
return Q;
}
-// ── Kalman predict step ──────────────────────────────────────────────────────
-void HarmonicGuider::kalmanPredict(double dt, double ra_pulse_px, double dec_pulse_px,
- double alt_deg, double parallactic_angle_deg)
+// ── Drift model (px/s); altitude clamped to the fitted range ─────────────────
+void HarmonicGuider::driftRates(double alt_deg, double parallactic_angle_deg,
+ double &ra_rate, double &dec_rate) const
+{
+ 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;
+
+ ra_rate = m_drift_ra;
+ dec_rate = m_drift_dec + m_d_polar;
+ if (std::abs(cos_alt) > 1e-4)
+ {
+ ra_rate += m_k_ref / (cos_alt * cos_alt);
+ dec_rate += m_k_ref_dec * std::sin(q_rad) / (cos_alt * cos_alt);
+ }
+}
+
+// ── Kalman predict step: free dynamics only, pulses are applied in update() ──
+void HarmonicGuider::kalmanPredict(double dt, double alt_deg, double parallactic_angle_deg)
{
- // Build state transition matrix
Eigen::Matrix<double, N_STATES, N_STATES> F;
buildF(F, dt);
- // Compute effective pulse (what actually moves the mount)
- // The spring absorbs κ fraction; only (1-κ) is immediately effective.
- // The absorbed part goes into the spring state.
- const double effective_ra = ra_pulse_px * (1.0 - m_kappa_ra);
- const double effective_dec = dec_pulse_px * (1.0 - m_kappa_dec);
-
- // Capture the spring tension BEFORE F decays it. The energy released to position
- // this frame is (spring_before - spring_after) = spring_before·(1 - e^(-dt/τ)).
- // Summed over all frames this telescopes to the full absorbed κ·pulse, so the
- // mount fully recovers the absorbed pulse — matching the training step response
- // pos(t) = pulse·(1 - κ·e^(-t/τ)) that train_harmonic.py fits κ/τ to (asymptote = pulse).
- // (The old code used the POST-decay spring value, which only ever released
- // κ·pulse·e^(-dt/τ) in total, permanently losing the rest.)
+ // Capture the spring tension BEFORE F decays it: released = spring_before·(1 - e^(-dt/τ)).
const double spring_before_ra = m_x(RA_SPRING);
const double spring_before_dec = m_x(DEC_SPRING);
- // Predict state (this decays the spring by e^(-dt/τ) via F)
m_x = F * m_x;
- // Add the spring absorption: new spring tension += κ * raw_pulse.
- // Freshly-absorbed energy releases on FUTURE frames, not this one, so it is added
- // after spring_before was captured.
- m_x(RA_SPRING) += m_kappa_ra * ra_pulse_px;
- m_x(DEC_SPRING) += m_kappa_dec * dec_pulse_px;
-
- // Add the effective pulse correction to position
- // (Negative because a correction pulse reduces the error)
- m_x(RA_POS) -= effective_ra;
- m_x(DEC_POS) -= effective_dec;
-
- // Release the spring into position. The release is part of the SAME correction as the
- // immediate term above (the mount continues moving in the correcting direction as the
- // spring unwinds), so it also SUBTRACTS from the error. Immediate (1-κ)·pulse plus the
- // total release κ·pulse then sum to the full pulse, matching the trainer step response
- // pos(t)=pulse·(1-κ·e^(-t/τ)) whose asymptote is the full pulse. released = spring_before - spring_after.
+ // Spring release continues the correction as absorbed pulse energy unwinds
const double spring_released_ra = spring_before_ra * (1.0 - std::exp(-dt / m_tau_ra));
const double spring_released_dec = spring_before_dec * (1.0 - std::exp(-dt / m_tau_dec));
m_x(RA_POS) -= spring_released_ra;
m_x(DEC_POS) -= spring_released_dec;
- // Add drift terms
- const double alt_rad = alt_deg * M_PI / 180.0;
- const double cos_alt = std::cos(alt_rad);
- const double q_rad = parallactic_angle_deg * M_PI / 180.0;
-
- double ra_drift_rate = m_drift_ra;
- double dec_drift_rate = m_drift_dec + m_d_polar;
- if (std::abs(cos_alt) > 1e-4)
- {
- ra_drift_rate += m_k_ref / (cos_alt * cos_alt);
- dec_drift_rate += m_k_ref_dec * std::sin(q_rad) / (cos_alt * cos_alt);
- }
-
+ double ra_drift_rate = 0.0, dec_drift_rate = 0.0;
+ driftRates(alt_deg, parallactic_angle_deg, ra_drift_rate, dec_drift_rate);
m_x(RA_POS) += ra_drift_rate * dt;
m_x(DEC_POS) += dec_drift_rate * dt;
@@ -405,7 +376,7 @@ void HarmonicGuider::kalmanPredict(double dt, double ra_pulse_px, double dec_pul
}
// ── Kalman update step ───────────────────────────────────────────────────────
-void HarmonicGuider::kalmanUpdate(double ra_meas_px, double dec_meas_px)
+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
@@ -418,10 +389,10 @@ void HarmonicGuider::kalmanUpdate(double ra_meas_px, double dec_meas_px)
H(1, DEC_PE_SIN) = 1.0;
}
- // Measurement noise (fixed from warmup SNR estimate)
- const double r_base = 0.1; // Base measurement noise in px²
- const double snr_factor = (m_lastSNR > 10.0) ? (100.0 / (m_lastSNR * m_lastSNR)) : 1.0;
- Eigen::Matrix2d R = Eigen::Matrix2d::Identity() * (r_base * snr_factor);
+ // Measurement noise from the current frame's SNR (~0.5 px at SNR 30)
+ const double snr_safe = std::max(snr, 5.0);
+ const double r_val = std::clamp(0.25 * (30.0 / snr_safe) * (30.0 / snr_safe), 0.04, 4.0);
+ Eigen::Matrix2d R = Eigen::Matrix2d::Identity() * r_val;
// Innovation
Eigen::Vector2d z(ra_meas_px, dec_meas_px);
@@ -465,18 +436,6 @@ GuideOutput HarmonicGuider::predict(const GuideFrameData &frame)
m_lastRAPulseMs = frame.ra_pulse_ms;
m_lastDECPulseMs = frame.dec_pulse_ms;
- // Convert pulse from ms to pixels for the Kalman filter:
- // pulse_px = pulse_ms / ms_per_arcsec / pixel_scale
- // ms_per_arcsec is the real guide-rate calibration forwarded from gmath — the SAME factor the
- // measurement path uses to remove this pulse. Fall back to a 1000 ms/arcsec placeholder only
- // when calibration isn't available yet (e.g. before the first calibration completes).
- const double ra_ms_per_arcsec = (frame.ra_ms_per_arcsec > 0.0) ? frame.ra_ms_per_arcsec : 1000.0;
- const double dec_ms_per_arcsec = (frame.dec_ms_per_arcsec > 0.0) ? frame.dec_ms_per_arcsec : 1000.0;
- const double ra_pulse_arcsec = std::abs(frame.ra_pulse_ms) / ra_ms_per_arcsec;
- const double dec_pulse_arcsec = std::abs(frame.dec_pulse_ms) / dec_ms_per_arcsec;
- const double ra_pulse_px = (frame.ra_pulse_ms >= 0 ? 1.0 : -1.0) * ra_pulse_arcsec / frame.pixel_scale;
- const double dec_pulse_px = (frame.dec_pulse_ms >= 0 ? 1.0 : -1.0) * dec_pulse_arcsec / frame.pixel_scale;
-
// Q-net input feature: raw frame-to-frame change of tracking error, computed the
// same way as train_harmonic.py (|ra_raw_px - prev_ra_raw_px|). Feeding the Kalman
// innovation here instead is out-of-distribution for the
@@ -494,19 +453,28 @@ GuideOutput HarmonicGuider::predict(const GuideFrameData &frame)
m_prevDecRawPx = frame.dec_raw_px;
m_hasPrevRaw = true;
- // Run Kalman predict step
- kalmanPredict(frame.dt, ra_pulse_px, dec_pulse_px,
- frame.altitude_deg, frame.parallactic_angle_deg);
+ // 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);
+ }
+
+ kalmanPredict(frame.dt, frame.altitude_deg, frame.parallactic_angle_deg);
- // The prediction for the next frame is the predicted position state
- // (which includes drift + PE + spring release effects)
- m_lastPredRA = m_x(RA_POS);
- m_lastPredDEC = m_x(DEC_POS);
+ double pred_ra = m_x(RA_POS);
+ double pred_dec = m_x(DEC_POS);
if (m_pe_period > 0.0)
{
- m_lastPredRA += m_x(RA_PE_SIN);
- m_lastPredDEC += m_x(DEC_PE_SIN);
+ pred_ra += m_x(RA_PE_SIN);
+ pred_dec += m_x(DEC_PE_SIN);
}
+
+ // Prediction is the expected uncorrected drift over the next interval
+ m_lastPredRA = pred_ra - post_ra;
+ m_lastPredDEC = pred_dec - post_dec;
m_hasLastPred = true;
GuideOutput out;
@@ -520,7 +488,7 @@ 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 = drift, mlp = spring + PE (the "learned" part)
+ // 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;
@@ -532,10 +500,20 @@ GuideOutput HarmonicGuider::predict(const GuideFrameData &frame)
}
void HarmonicGuider::update(double /*ra_error_px*/, double /*dec_error_px*/,
- double uncorrected_drift_ra_px, double uncorrected_drift_dec_px, double snr)
+ double uncorrected_drift_ra_px, double uncorrected_drift_dec_px, double snr,
+ double ra_pulse_px, double dec_pulse_px)
{
- // Run Kalman update with actual measurement
- kalmanUpdate(uncorrected_drift_ra_px, uncorrected_drift_dec_px);
+ // The absorbed κ·pulse enters the uncorrected trajectory now and releases later
+ m_x(RA_SPRING) += m_kappa_ra * ra_pulse_px;
+ m_x(DEC_SPRING) += m_kappa_dec * dec_pulse_px;
+ m_x(RA_POS) += m_kappa_ra * ra_pulse_px;
+ m_x(DEC_POS) += m_kappa_dec * dec_pulse_px;
+
+ // The filter observes the integrated uncorrected position
+ m_uncorrPosRA += uncorrected_drift_ra_px;
+ m_uncorrPosDEC += uncorrected_drift_dec_px;
+
+ kalmanUpdate(m_uncorrPosRA, m_uncorrPosDEC, snr);
if (m_hasLastPred)
{
@@ -600,34 +578,53 @@ void HarmonicGuider::updateConfidence(double innovRA, double innovDec, double sn
m_confidence = std::clamp(warmup_factor * snr_factor * prediction_quality, 0.0, 1.0);
}
-// ── Dark guiding prediction ──────────────────────────────────────────────────
+// ── Dark guiding: propagate a copy statelessly, return the interval increment ─
GuideOutput HarmonicGuider::darkPredict(double dt_sec)
{
m_lastSessionSec += dt_sec;
- // Run Kalman predict without any pulse input
- kalmanPredict(dt_sec, 0.0, 0.0,
- m_lastAltRad * 180.0 / M_PI, m_lastParallacticAngleDeg);
+ Eigen::Matrix<double, N_STATES, 1> x = m_x;
+ Eigen::Matrix<double, N_STATES, N_STATES> F;
+ buildF(F, dt_sec);
- double pred_ra = m_x(RA_POS);
- double pred_dec = m_x(DEC_POS);
+ double post_ra = x(RA_POS);
+ double post_dec = x(DEC_POS);
if (m_pe_period > 0.0)
{
- pred_ra += m_x(RA_PE_SIN);
- pred_dec += m_x(DEC_PE_SIN);
+ post_ra += x(RA_PE_SIN);
+ post_dec += x(DEC_PE_SIN);
+ }
+
+ 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));
+ x = F * x;
+ x(RA_POS) -= release_ra;
+ x(DEC_POS) -= release_dec;
+
+ double ra_rate = 0.0, dec_rate = 0.0;
+ driftRates(m_lastAltRad * 180.0 / M_PI, m_lastParallacticAngleDeg, ra_rate, dec_rate);
+ 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);
}
GuideOutput out;
out.valid = (m_frameCount > warmupFrames());
out.confidence = m_confidence;
- out.ra_correction_arcsec = pred_ra * m_lastPixelScale;
- out.dec_correction_arcsec = pred_dec * m_lastPixelScale;
+ out.ra_correction_arcsec = (pred_ra - post_ra) * m_lastPixelScale;
+ out.dec_correction_arcsec = (pred_dec - post_dec) * m_lastPixelScale;
- out.physics_ra_arcsec = m_x(RA_VEL) * dt_sec * m_lastPixelScale;
- out.physics_dec_arcsec = m_x(DEC_VEL) * dt_sec * m_lastPixelScale;
- out.mlp_ra_arcsec = (pred_ra - m_x(RA_VEL) * dt_sec) * m_lastPixelScale;
- out.mlp_dec_arcsec = (pred_dec - m_x(DEC_VEL) * dt_sec) * m_lastPixelScale;
+ 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;
return out;
}
diff --git a/kstars/ekos/guide/internalguide/harmonic_guider.h b/kstars/ekos/guide/internalguide/harmonic_guider.h
index 40cd26356b..49b4b16dd8 100644
--- a/kstars/ekos/guide/internalguide/harmonic_guider.h
+++ b/kstars/ekos/guide/internalguide/harmonic_guider.h
@@ -34,7 +34,8 @@ class HarmonicGuider : public MountSpecificGuider
GuideOutput predict(const GuideFrameData &frame) override;
GuideOutput darkPredict(double dt_sec) override;
void update(double ra_error_px, double dec_error_px, double uncorrected_drift_ra_px,
- double uncorrected_drift_dec_px, double snr) override;
+ double uncorrected_drift_dec_px, double snr,
+ double ra_pulse_px, double dec_pulse_px) override;
double confidence() const override
{
return m_confidence;
@@ -82,10 +83,14 @@ class HarmonicGuider : public MountSpecificGuider
double m_k_ref { 0.0 }; ///< RA refraction coefficient
double m_d_polar { 0.0 }; ///< DEC polar drift rate (px/s)
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 };
// ── Kalman filter state (static to survive object recreation) ────────
static Eigen::Matrix<double, N_STATES, 1> m_x; ///< State estimate
static Eigen::Matrix<double, N_STATES, N_STATES> m_P; ///< Error covariance
+ static double m_uncorrPosRA; ///< Integrated uncorrected position (the measurement)
+ 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
@@ -137,8 +142,9 @@ class HarmonicGuider : public MountSpecificGuider
void buildF(Eigen::Matrix<double, N_STATES, N_STATES> &F, double dt) const;
Eigen::Matrix<double, N_STATES, N_STATES> computeQ(double snr, double snr_delta,
double innov_ra, double innov_dec, double dt) const;
- void kalmanPredict(double dt, double ra_pulse_px, double dec_pulse_px,
- double alt_deg, double parallactic_angle_deg);
- void kalmanUpdate(double ra_meas_px, double dec_meas_px);
+ void driftRates(double alt_deg, double parallactic_angle_deg,
+ double &ra_rate, double &dec_rate) const;
+ 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);
};
diff --git a/kstars/ekos/guide/internalguide/mount_guider.h b/kstars/ekos/guide/internalguide/mount_guider.h
index fe7d483e71..7b75ce1a7e 100644
--- a/kstars/ekos/guide/internalguide/mount_guider.h
+++ b/kstars/ekos/guide/internalguide/mount_guider.h
@@ -106,9 +106,12 @@ class MountSpecificGuider
* @param uncorrected_drift_ra_px Uncorrected physical RA drift (pixels)
* @param uncorrected_drift_dec_px Uncorrected physical DEC drift (pixels)
* @param snr Guide star SNR
+ * @param ra_pulse_px Signed RA pulse applied over this interval (pixels)
+ * @param dec_pulse_px Signed DEC pulse applied over this interval (pixels)
*/
virtual void update(double ra_error_px, double dec_error_px, double uncorrected_drift_ra_px,
- double uncorrected_drift_dec_px, double snr) = 0;
+ double uncorrected_drift_dec_px, double snr,
+ double ra_pulse_px, double dec_pulse_px) = 0;
/**
* @brief Current model confidence in [0, 1].
diff --git a/kstars/ekos/guide/internalguide/worm_gear_guider.cpp b/kstars/ekos/guide/internalguide/worm_gear_guider.cpp
index 6b2345302c..e7b64af198 100644
--- a/kstars/ekos/guide/internalguide/worm_gear_guider.cpp
+++ b/kstars/ekos/guide/internalguide/worm_gear_guider.cpp
@@ -7,6 +7,7 @@
#include "worm_gear_guider.h"
#include "Options.h"
+#include "ekos_guide_debug.h"
#include <QFile>
#include <QJsonDocument>
@@ -47,60 +48,47 @@ bool WormGearGuider::validateFingerprint(const QJsonObject &fp)
if (fp.isEmpty())
return true;
- if (fp.contains("guide_exposure_s"))
+ const struct
{
- const double expected = fp["guide_exposure_s"].toDouble();
- if (!fpDoubleClose(expected, Options::guideExposure(), 0.05))
- return false;
- }
-
- if (fp.contains("ra_proportional_gain") &&
- !fpDoubleClose(fp["ra_proportional_gain"].toDouble(), Options::rAProportionalGain()))
- return false;
-
- if (fp.contains("dec_proportional_gain") &&
- !fpDoubleClose(fp["dec_proportional_gain"].toDouble(), Options::dECProportionalGain()))
- return false;
-
- if (fp.contains("ra_integral_gain") &&
- !fpDoubleClose(fp["ra_integral_gain"].toDouble(), Options::rAIntegralGain()))
- return false;
-
- if (fp.contains("dec_integral_gain") &&
- !fpDoubleClose(fp["dec_integral_gain"].toDouble(), Options::dECIntegralGain()))
- return false;
-
- if (fp.contains("ra_min_pulse_arcsec") &&
- !fpDoubleClose(fp["ra_min_pulse_arcsec"].toDouble(), Options::rAMinimumPulseArcSec()))
- return false;
-
- if (fp.contains("dec_min_pulse_arcsec") &&
- !fpDoubleClose(fp["dec_min_pulse_arcsec"].toDouble(), Options::dECMinimumPulseArcSec()))
- return false;
-
- if (fp.contains("ra_max_pulse_arcsec") &&
- !fpDoubleClose(fp["ra_max_pulse_arcsec"].toDouble(), Options::rAMaximumPulseArcSec()))
- return false;
-
- if (fp.contains("dec_max_pulse_arcsec") &&
- !fpDoubleClose(fp["dec_max_pulse_arcsec"].toDouble(), Options::dECMaximumPulseArcSec()))
- return false;
-
- if (fp.contains("ra_hysteresis") &&
- !fpDoubleClose(fp["ra_hysteresis"].toDouble(), Options::rAHysteresis()))
- return false;
+ const char *key;
+ double current;
+ double tol;
+ } checks[] =
+ {
+ { "guide_exposure_s", Options::guideExposure(), 0.05 },
+ { "ra_proportional_gain", Options::rAProportionalGain(), 1e-4 },
+ { "dec_proportional_gain", Options::dECProportionalGain(), 1e-4 },
+ { "ra_integral_gain", Options::rAIntegralGain(), 1e-4 },
+ { "dec_integral_gain", Options::dECIntegralGain(), 1e-4 },
+ { "ra_min_pulse_arcsec", Options::rAMinimumPulseArcSec(), 1e-4 },
+ { "dec_min_pulse_arcsec", Options::dECMinimumPulseArcSec(), 1e-4 },
+ { "ra_max_pulse_arcsec", static_cast<double>(Options::rAMaximumPulseArcSec()), 1e-4 },
+ { "dec_max_pulse_arcsec", static_cast<double>(Options::dECMaximumPulseArcSec()), 1e-4 },
+ { "ra_hysteresis", Options::rAHysteresis(), 1e-4 },
+ { "dec_hysteresis", Options::dECHysteresis(), 1e-4 },
+ };
- if (fp.contains("dec_hysteresis") &&
- !fpDoubleClose(fp["dec_hysteresis"].toDouble(), Options::dECHysteresis()))
- return false;
+ bool ok = true;
+ for (const auto &c : checks)
+ {
+ if (fp.contains(c.key) && !fpDoubleClose(fp[c.key].toDouble(), c.current, c.tol))
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected:" << c.key << "recorded"
+ << fp[c.key].toDouble() << "current" << c.current;
+ ok = false;
+ }
+ }
- if (fp.contains("guide_binning") &&
- fp["guide_binning"].toString() != Options::guideBinning())
- return false;
+ if (fp.contains("guide_binning") && fp["guide_binning"].toString() != Options::guideBinning())
+ {
+ qCWarning(KSTARS_EKOS_GUIDE) << "AI weights rejected: guide_binning recorded"
+ << fp["guide_binning"].toString() << "current" << Options::guideBinning();
+ ok = false;
+ }
// Pulse algorithm is Standard (0) in training fingerprint; runtime uses AI — skip.
- return true;
+ return ok;
}
bool WormGearGuider::loadWeights(const QString &weightsPath)
@@ -130,6 +118,8 @@ bool WormGearGuider::loadWeights(const QString &weightsPath)
m_d_ra_extra = phys["d_ra_extra"].toDouble(0.0);
m_d_polar = phys["d_polar"].toDouble(0.0);
m_k_ref_dec = phys["k_ref_dec"].toDouble(0.0);
+ m_fit_alt_min = phys["fit_alt_min"].toDouble(35.0);
+ m_fit_alt_max = phys["fit_alt_max"].toDouble(65.0);
QJsonObject norm = root["normalization"].toObject();
m_alt_scale = norm["alt_scale"].toDouble(90.0);
@@ -266,7 +256,8 @@ GuideOutput WormGearGuider::predict(const GuideFrameData &frame)
}
void WormGearGuider::update(double /*ra_error_px*/, double /*dec_error_px*/,
- double uncorrected_drift_ra_px_delta, double uncorrected_drift_dec_px_delta, double snr)
+ double uncorrected_drift_ra_px_delta, double uncorrected_drift_dec_px_delta, double snr,
+ double /*ra_pulse_px*/, double /*dec_pulse_px*/)
{
m_uncorrectedPosRA += uncorrected_drift_ra_px_delta;
m_uncorrectedPosDEC += uncorrected_drift_dec_px_delta;
@@ -304,7 +295,8 @@ double WormGearGuider::physicsRA(double t_sec, double /*altitude_deg*/) const
double WormGearGuider::physicsDEC(double altitude_deg, double parallactic_angle_deg) const
{
- const double alt_rad = altitude_deg * M_PI / 180.0;
+ // Refraction fit is only valid inside the fitted altitude range.
+ const double alt_rad = std::clamp(altitude_deg, m_fit_alt_min, m_fit_alt_max) * M_PI / 180.0;
const double q_rad = parallactic_angle_deg * M_PI / 180.0;
const double cos_alt = std::cos(alt_rad);
diff --git a/kstars/ekos/guide/internalguide/worm_gear_guider.h b/kstars/ekos/guide/internalguide/worm_gear_guider.h
index f6820309e6..341cf0bfb5 100644
--- a/kstars/ekos/guide/internalguide/worm_gear_guider.h
+++ b/kstars/ekos/guide/internalguide/worm_gear_guider.h
@@ -33,7 +33,7 @@ class WormGearGuider : public MountSpecificGuider
GuideOutput predict(const GuideFrameData &frame) override;
GuideOutput darkPredict(double dt_sec) override;
void update(double ra_error_px, double dec_error_px, double uncorrected_drift_ra_px, double uncorrected_drift_dec_px,
- double snr) override;
+ double snr, double ra_pulse_px, double dec_pulse_px) override;
double confidence() const override
{
return m_confidence;
@@ -56,6 +56,8 @@ class WormGearGuider : public MountSpecificGuider
static double m_d_ra_extra; ///< Continuous RA drift rate (pixels/second)
double m_d_polar { 0.0 }; ///< Polar drift rate (pixels/second)
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 };
// ── Online phase estimation (4-State Position RLS) ────────────────
static Eigen::Vector4d m_rls_theta; ///< [sin_coeff, cos_coeff, v, C]
diff --git a/kstars/ekos/guide/offlinetrainer/train_direct_drive.py b/kstars/ekos/guide/offlinetrainer/train_direct_drive.py
index 1a22f096ba..dc2a6a9d08 100644
--- a/kstars/ekos/guide/offlinetrainer/train_direct_drive.py
+++ b/kstars/ekos/guide/offlinetrainer/train_direct_drive.py
@@ -142,6 +142,8 @@ def train_direct_drive(sysid: dict,
"k_ref_dec": float(k_ref_dec),
"d_ra_extra": float(d_ra_extra),
"phi_drift": float(phi_drift),
+ "fit_alt_min": float(np.min(altitudes)),
+ "fit_alt_max": float(np.max(altitudes)),
},
"training_stats": {
"n_measurements": int(len(ra_drifts)),
diff --git a/kstars/ekos/guide/offlinetrainer/train_harmonic.py b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
index 85bae93e8b..2982a3dd2b 100644
--- a/kstars/ekos/guide/offlinetrainer/train_harmonic.py
+++ b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
@@ -66,6 +66,12 @@ def train_harmonic(sysid: dict,
drift_ra, drift_dec, d_polar, k_ref, k_ref_dec = _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]
+ fit_alt_min = min(fit_alts) if fit_alts else 35.0
+ fit_alt_max = max(fit_alts) if fit_alts else 65.0
+
if verbose:
print(f" drift_ra={drift_ra:.6e} px/s drift_dec={drift_dec:.6e} px/s")
print(f" d_polar={d_polar:.6e} px/s")
@@ -96,6 +102,8 @@ def train_harmonic(sysid: dict,
"k_ref": float(k_ref),
"d_polar": float(d_polar),
"k_ref_dec": float(k_ref_dec),
+ "fit_alt_min": float(fit_alt_min),
+ "fit_alt_max": float(fit_alt_max),
},
"qnet": qnet_weights,
}
@@ -105,26 +113,17 @@ def train_harmonic(sysid: dict,
# Phase 1: Spring parameter fitting from pulse_response sessions
# ═══════════════════════════════════════════════════════════════════════════════
-def _spring_response_model(t, pulse_magnitude, kappa, tau):
- """
- Model: the mount position after a pulse = pulse_mag * (1 - κ * exp(-t/τ))
- At t=0 the immediate response is pulse_mag * (1 - κ).
- As t→∞ the spring releases and total response → pulse_mag.
- """
- return pulse_magnitude * (1.0 - kappa * np.exp(-t / tau))
-
-
def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
"""
Fit spring constant κ and time constant τ from pulse_response sessions.
- For each pulse_response session matching the given axis:
- 1. Extract the response curve (position change after the pulse)
- 2. Fit: response(t) = pulse_mag * (1 - κ * exp(-t/τ))
- 3. Average κ and τ across all pulse magnitudes for robustness
+ 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)
"""
+ # Unmeasured means unmodeled: the default kappa stays 0
+ DEFAULTS = (0.0, 1.5)
pulse_sessions = [
s for s in sysid["sessions"]
if s.get("type") == "pulse_response" and s.get("pulse_axis", "").upper() == axis.upper()
@@ -133,76 +132,171 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
if not pulse_sessions:
if verbose:
print(f" [{axis}] No pulse_response sessions found. Using defaults (κ=0.2, τ=1.5s)")
- return 0.2, 1.5
+ return DEFAULTS
- kappas = []
- taus = []
+ axis_key = "ra_raw_px" if axis.upper() == "RA" else "dec_raw_px"
- for s in pulse_sessions:
- pulse_mag = s.get("pulse_magnitude_ms", 100.0)
+ def session_curve(s):
+ """(t, signed displacement from baseline) for one pulse session, or None."""
frames = s.get("response_frames", [])
- if len(frames) < 3:
- continue
-
- # Build time and position arrays
- # Position is the cumulative displacement from the first frame
- axis_key = "ra_raw_px" if axis.upper() == "RA" else "dec_raw_px"
- baseline = frames[0].get(axis_key, 0.0)
-
- t_vals = []
- pos_vals = []
- for i, f in enumerate(frames):
- if i == 0:
- continue # Skip the frame where pulse was sent
- t = f.get("t", 0.0) - frames[0].get("t", 0.0)
- if t <= 0:
- t = i * guide_exp
- t_vals.append(t)
- pos_vals.append(abs(f.get(axis_key, 0.0) - baseline))
-
- t_arr = np.array(t_vals)
- pos_arr = np.array(pos_vals)
-
- if len(t_arr) < 3 or np.max(pos_arr) < 0.01:
- continue
-
- # Normalize position by the expected full response
- # (we don't know the exact calibration, so use max observed position)
- max_response = np.max(pos_arr)
-
+ 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:
- # Fit the spring model
- # response(t) = max_response * (1 - κ * exp(-t/τ))
- def model(t, kappa, tau):
- return max_response * (1.0 - kappa * np.exp(-t / tau))
-
- popt, pcov = scipy.optimize.curve_fit(
- model, t_arr, pos_arr,
- p0=[0.3, 1.5],
- bounds=([0.0, 0.1], [0.9, 10.0]),
- maxfev=5000
- )
- kappas.append(popt[0])
- taus.append(popt[1])
+ 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
+
+ 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])
if verbose:
- print(f" [{axis}] Pulse {pulse_mag}ms "
- f"{s.get('pulse_direction', '?')}: "
- f"κ={popt[0]:.3f}, τ={popt[1]:.2f}s "
- f"(max_response={max_response:.3f}px)")
- except (RuntimeError, ValueError) as e:
+ 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)))
if verbose:
- print(f" [{axis}] Pulse {pulse_mag}ms: curve_fit failed ({e})")
- continue
+ 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}] All curve fits failed. Using defaults.")
- return 0.2, 1.5
+ 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 DEFAULTS
+
+ # 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 DEFAULTS
- # Median is more robust than mean against outliers
kappa_result = float(np.median(kappas))
- tau_result = float(np.median(taus))
+ tau_result = float(np.median(taus)) if taus 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 DEFAULTS
if verbose:
print(f" [{axis}] Final: κ={kappa_result:.3f} (from {len(kappas)} fits), "
@@ -215,86 +309,118 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
# Phase 2: PE period detection from free-drift data
# ═══════════════════════════════════════════════════════════════════════════════
+def _pulse_correction_px(pulse_ms, cal_rate_ms_per_arcsec, pixel_scale):
+ """Signed pulse displacement in pixels (same helper as train_worm_gear)."""
+ if pulse_ms == 0.0 or cal_rate_ms_per_arcsec <= 0.0:
+ return 0.0
+ return (pulse_ms / cal_rate_ms_per_arcsec) / pixel_scale
+
+
+def _pe_candidate_series(sysid: dict, guide_exp: float):
+ """
+ Build (t, position, label) series for PE search.
+
+ Free drifts give the uncorrected trajectory directly. Standard-guiding
+ sessions are much longer (8 min vs 2 min) and reach the wave-generator
+ periods (~300-900s), so reconstruct their uncorrected trajectory by
+ adding the applied pulses back (same compensated formula as train_worm_gear:
+ JSON RA pulses are ADDED, DEC pulses SUBTRACTED).
+ """
+ pixel_scale = sysid["equipment"].get("pixel_scale_arcsec_per_px", 1.0)
+ series = []
+ for s in sysid["sessions"]:
+ frames = s.get("frames", [])
+ if len(frames) < 20:
+ continue
+ if s["type"] == "free_drift":
+ t, pos = 0.0, []
+ t_vals = []
+ for f in frames:
+ t += f.get("dt", guide_exp)
+ t_vals.append(t)
+ pos.append(f["ra_raw_px"])
+ series.append((np.array(t_vals), np.array(pos), f"free_drift {s.get('session_id', '?')}"))
+ elif s["type"] == "standard_guiding":
+ cal = s.get("ra_ms_per_arcsec", 0.0)
+ if cal <= 0.0 or pixel_scale <= 0.0:
+ continue
+ t, p = 0.0, 0.0
+ t_vals, pos = [], []
+ for i in range(1, len(frames)):
+ t += frames[i].get("dt", guide_exp)
+ p += (frames[i]["ra_raw_px"] - frames[i - 1]["ra_raw_px"]
+ + _pulse_correction_px(frames[i - 1].get("ra_pulse_ms", 0.0), cal, pixel_scale))
+ t_vals.append(t)
+ pos.append(p)
+ series.append((np.array(t_vals), np.array(pos), f"standard {s.get('session_id', '?')}"))
+ return series
+
+
def _estimate_pe(sysid: dict, guide_exp: float, verbose: bool):
"""
- Detect PE period and amplitude from free-drift data using Lomb-Scargle.
- Searches 0.1-0.5 Hz (periods 2-10s) for harmonic drive PE.
- Returns (period_seconds, amplitude_pixels) or (0.0, 0.0) if no PE detected.
+ Detect PE period and amplitude via Lomb-Scargle over every usable session.
+ Band per series: from 2 cycles per observation span up to Nyquist, so long
+ standard-guiding sessions expose the strain-wave fundamental (~300-900s)
+ and short free drifts still cover fast components.
+ Returns (period_seconds, amplitude_pixels) or (0.0, 0.0) if none significant.
"""
- free_drift_sessions = [s for s in sysid["sessions"] if s["type"] == "free_drift"]
- if not free_drift_sessions:
+ candidates = _pe_candidate_series(sysid, guide_exp)
+ if not candidates:
if verbose:
- print(" No free_drift sessions found. PE detection skipped.")
+ print(" No usable sessions found. PE detection skipped.")
return 0.0, 0.0
- # Use the longest free drift session
- longest_session = max(free_drift_sessions, key=lambda s: len(s["frames"]))
- frames = longest_session["frames"]
+ best = None # (snr, period, amplitude, label)
+ for t_arr, ra_arr, label in candidates:
+ span = t_arr[-1] - t_arr[0]
+ nyquist = 0.5 / guide_exp
+ f_min = max(2.0 / span, 0.002) # need >= 2 observed cycles
+ f_max = min(0.5, nyquist * 0.9)
+ if f_min >= f_max or span <= 0:
+ continue
- t_vals = []
- ra_vals = []
- t = 0.0
- for f in frames:
- t += f.get("dt", guide_exp)
- t_vals.append(t)
- ra_vals.append(f["ra_raw_px"])
+ slope, intercept, _, _, _ = scipy.stats.linregress(t_arr, ra_arr)
+ ra_detrended = ra_arr - (slope * t_arr + intercept)
- t_arr = np.array(t_vals)
- ra_arr = np.array(ra_vals)
+ f_search = np.geomspace(f_min, f_max, 4000)
+ omega = 2 * np.pi * f_search
+ Pxx = scipy.signal.lombscargle(t_arr, ra_detrended, omega, precenter=True)
- if len(t_arr) < 20:
- if verbose:
- print(" Free drift session too short for PE detection.")
- return 0.0, 0.0
+ peak_idx = np.argmax(Pxx)
+ peak_freq = f_search[peak_idx]
+ noise_floor = np.median(Pxx)
+ snr = Pxx[peak_idx] / (noise_floor + 1e-10)
+ amplitude = np.sqrt(4 * Pxx[peak_idx] / len(t_arr))
- # Detrend to remove linear drift
- slope, intercept, _, _, _ = scipy.stats.linregress(t_arr, ra_arr)
- ra_detrended = ra_arr - (slope * t_arr + intercept)
+ at_edge = peak_freq <= f_min * 1.05
+ if verbose:
+ print(f" [LS] {label}: span={span:.0f}s, band {1/f_max:.1f}-{1/f_min:.0f}s, "
+ f"peak {1/peak_freq:.1f}s, SNR {snr:.1f}, amp {amplitude:.3f}px"
+ f"{' [AT BAND EDGE]' if at_edge else ''}")
- # Lomb-Scargle in the harmonic drive PE frequency range
- # PE periods 2-10s → frequencies 0.1-0.5 Hz
- nyquist = 0.5 / guide_exp # Maximum observable frequency
- f_max = min(0.5, nyquist * 0.9) # Stay below Nyquist
- f_min = 0.1
+ if best is None or snr > best[0]:
+ best = (snr, 1.0 / peak_freq, amplitude, label, at_edge)
- if f_min >= f_max:
+ if best is None or best[0] < 10.0:
if verbose:
- print(f" Guide exposure ({guide_exp}s) too long for harmonic PE detection. "
- f"Nyquist={nyquist:.2f} Hz")
+ snr = 0.0 if best is None else best[0]
+ print(f" [LS] PE not significant (best SNR {snr:.1f} < 10.0). Disabling PE states.")
return 0.0, 0.0
- f_search = np.linspace(f_min, f_max, 2000)
- omega = 2 * np.pi * f_search
- Pxx = scipy.signal.lombscargle(t_arr, ra_detrended, omega, precenter=True)
-
- peak_idx = np.argmax(Pxx)
- peak_power = Pxx[peak_idx]
- peak_freq = f_search[peak_idx]
- peak_period = 1.0 / peak_freq
-
- # Significance test: is the peak significantly above the noise floor?
- # Use the median power as the noise level estimate
- noise_floor = np.median(Pxx)
- snr = peak_power / (noise_floor + 1e-10)
-
- # Amplitude estimate from Lomb-Scargle power
- amplitude = np.sqrt(4 * peak_power / len(t_arr))
-
- if verbose:
- print(f" [LS] Searched {f_min:.2f}-{f_max:.2f} Hz ({1/f_max:.1f}-{1/f_min:.1f}s)")
- print(f" [LS] Peak at {peak_freq:.4f} Hz → Period: {peak_period:.2f}s")
- print(f" [LS] Peak SNR: {snr:.1f} (threshold: 10.0)")
- print(f" [LS] Amplitude: {amplitude:.4f} px")
-
- # Require SNR > 10 for a significant PE detection
- if snr < 10.0:
+ snr, peak_period, amplitude, label, at_edge = best
+ if at_edge:
+ # True period is longer than the session can resolve: disable PE, tell the user
if verbose:
- print(f" [LS] PE not significant (SNR {snr:.1f} < 10.0). Disabling PE states.")
+ print(f" [LS] WARNING: dominant PE ({amplitude:.2f}px, SNR {snr:.0f}) sits at the "
+ f"band edge ({peak_period:.0f}s) — true period is longer than the session "
+ f"can resolve. Collect a session of at least {2.5 * peak_period:.0f}s "
+ f"(free drift or standard guiding) and retrain. Disabling PE states.")
return 0.0, 0.0
- # Sanity bounds
- peak_period = np.clip(peak_period, 1.5, 15.0)
- amplitude = np.clip(amplitude, 0.01, 5.0)
+ peak_period = np.clip(peak_period, 1.5, 1500.0)
+ amplitude = np.clip(amplitude, 0.01, 50.0)
+ if verbose:
+ print(f" [LS] Selected: {peak_period:.1f}s, amp {amplitude:.3f}px (from {label})")
return float(peak_period), float(amplitude)
diff --git a/kstars/ekos/guide/offlinetrainer/train_worm_gear.py b/kstars/ekos/guide/offlinetrainer/train_worm_gear.py
index dd5fd08450..8e05706b33 100644
--- a/kstars/ekos/guide/offlinetrainer/train_worm_gear.py
+++ b/kstars/ekos/guide/offlinetrainer/train_worm_gear.py
@@ -97,7 +97,11 @@ def train_worm_gear(sysid: dict,
"k_ref": float(k_ref),
"d_ra_extra": float(d_ra_extra),
"d_polar": float(d_polar),
- "k_ref_dec": float(k_ref_dec)
+ "k_ref_dec": float(k_ref_dec),
+ "fit_alt_min": float(min((s.get("altitude_deg", 45.0) for s in sysid["sessions"]
+ if s.get("type") == "free_drift"), default=35.0)),
+ "fit_alt_max": float(max((s.get("altitude_deg", 45.0) for s in sysid["sessions"]
+ if s.get("type") == "free_drift"), default=65.0))
},
"normalization": {
"alt_scale": 90.0,