Kstars

gpg.cpp
1/*
2 SPDX-FileCopyrightText: 2020 Hy Murveit <hy@murveit.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5*/
6
7#include "gpg.h"
8#include <vector>
9
10#include "guidestars.h"
11#include "Options.h"
12#include "ekos_guide_debug.h"
13
14#include "calibration.h"
15
16#include <QElapsedTimer>
17
18namespace
19{
20
21// Don't have GPG model large errors. Instantaneous large errors are probably not due to
22// periodic error, but rather one-offs.
23constexpr double MAX_ARCSEC_ERROR = 5.0;
24
25// If we skip this many samples, reset gpg.
26constexpr int MAX_SKIPPED_SAMPLES = 4;
27
28// Fills parameters from the KStars options.
29void getGPGParameters(GaussianProcessGuider::guide_parameters *parameters)
30{
31 // Parameters from standard options.
32 parameters->control_gain_ = Options::rAProportionalGain();
33 // We need a calibration to determine min move. This is re-set in computePulse() below.
34 parameters->min_move_ = 0.25;
35
36 // Parameters from GPG-specific options.
37 parameters->min_periods_for_inference_ = Options::gPGMinPeriodsForInference();
38 parameters->SE0KLengthScale_ = Options::gPGSE0KLengthScale();
39 parameters->SE0KSignalVariance_ = Options::gPGSE0KSignalVariance();
40 parameters->PKLengthScale_ = Options::gPGPKLengthScale();
41 parameters->PKPeriodLength_ = Options::gPGPeriod();
42 parameters->PKSignalVariance_ = Options::gPGPKSignalVariance();
43 parameters->SE1KLengthScale_ = Options::gPGSE1KLengthScale();
44 parameters->SE1KSignalVariance_ = Options::gPGSE1KSignalVariance();
45 parameters->min_periods_for_period_estimation_ = Options::gPGMinPeriodsForPeriodEstimate();
46 parameters->points_for_approximation_ = Options::gPGPointsForApproximation();
47 parameters->prediction_gain_ = Options::gPGpWeight();
48 parameters->compute_period_ = Options::gPGEstimatePeriod();
49}
50
51// Returns the SNR returned by guideStars, or if guideStars is null (e.g. we aren't
52// running SEP MultiStar, then just 50, as this isn't a critical parameter.
53double getSNR(const GuideStars *guideStars, const double raArcsecError)
54{
55 if (fabs(raArcsecError) > MAX_ARCSEC_ERROR)
56 return 1.0;
57
58 // If no SEP MultiStar, just fudge the SNR. Won't be a big deal.
59 if (guideStars == nullptr)
60 return 50.0;
61 else
62 return guideStars->getGuideStarSNR();
63}
64
65const QString directionStr(GuideDirection dir)
66{
67 switch (dir)
68 {
69 case RA_DEC_DIR:
70 return "Decrease RA";
71 case RA_INC_DIR:
72 return "Increase RA";
73 case DEC_DEC_DIR:
74 return "Decrease DEC";
75 case DEC_INC_DIR:
76 return "Increase DEC";
77 default:
78 return "NO DIR";
79 }
80}
81
82} // namespace
83
84double GPG::predictionContribution()
85{
86 if (gpg.get() == nullptr)
87 {
88 qCDebug(KSTARS_EKOS_GUIDE) << "Nullptr in predictionContribution()";
89 return 0.0;
90 }
91 return gpg->getPredictionContribution();
92}
93
94void GPG::updateParameters()
95{
96 // Parameters would be set when the gpg stars up.
97 if (gpg.get() == nullptr) return;
98
100 getGPGParameters(&parameters);
101 gpg->SetControlGain(parameters.control_gain_);
102 gpg->SetPeriodLengthsInference(parameters.min_periods_for_inference_);
103 gpg->SetMinMove(parameters.min_move_);
104 gpg->SetPeriodLengthsPeriodEstimation(parameters.min_periods_for_period_estimation_);
105 gpg->SetNumPointsForApproximation(parameters.points_for_approximation_);
106 gpg->SetPredictionGain(parameters.prediction_gain_);
107 gpg->SetBoolComputePeriod(parameters.compute_period_);
108
109 // The GPG header really should be in a namespace so NumParameters
110 // is not in a global namespace.
111 std::vector<double> hyperparameters(NumParameters);
112 hyperparameters[SE0KLengthScale] = parameters.SE0KLengthScale_;
113 hyperparameters[SE0KSignalVariance] = parameters.SE0KSignalVariance_;
114 hyperparameters[PKLengthScale] = parameters.PKLengthScale_;
115 hyperparameters[PKSignalVariance] = parameters.PKSignalVariance_;
116 hyperparameters[SE1KLengthScale] = parameters.SE1KLengthScale_;
117 hyperparameters[SE1KSignalVariance] = parameters.SE1KSignalVariance_;
118 hyperparameters[PKPeriodLength] = parameters.PKPeriodLength_;
119 gpg->SetGPHyperparameters(hyperparameters);
120}
121
122
123GPG::GPG()
124{
126 getGPGParameters(&parameters);
127 gpg.reset(new GaussianProcessGuider(parameters));
128 reset();
129}
130
131void GPG::reset()
132{
133 gpgSamples = 0;
134 gpgSkippedSamples = 0;
135 gpg->reset();
136 qCDebug(KSTARS_EKOS_GUIDE) << "Resetting GPG";
137}
138
139void GPG::startDithering(double dx, double dy, const Calibration &cal)
140{
141 // convert the x and y offsets to RA and DEC offsets.
142 double raPixels, decPixels;
143 cal.rotateToRaDec(dx, dy, &raPixels, &decPixels);
144
145 // The 1000 in the denominator below converts gear-milliseconds into gear-seconds,
146 // which is the unit GPG wants. Since calibration uses arcseconds, we convert using
147 // arcseconds/pixel, though that's inexact when the pixels aren't square.
148 // amount = (pixels * ms/pixel) / 1000;
149 const double amount = raPixels * cal.raPulseMillisecondsPerArcsecond() * cal.xArcsecondsPerPixel() / 1000.0;
150
151 qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither started. Gear-seconds =" << amount << "x,y was: " << dx << dy;
152 gpg->GuidingDithered(amount, 1.0);
153}
154
155void GPG::ditheringSettled(bool success)
156{
157 gpg->GuidingDitherSettleDone(success);
158 if (success)
159 qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither done (success)";
160 else
161 qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither done (failed)";
162}
163
164void GPG::suspended(const GuiderUtils::Vector &guideStarPosition,
165 const GuiderUtils::Vector &reticlePosition,
166 GuideStars *guideStars,
167 const Calibration &cal)
168{
169 constexpr int MaxGpgSamplesForReset = 25;
170 // We just reset the gpg if there's not enough samples to make
171 // it worth our while to keep it going. We're probably worse off
172 // to start it and then feed it bad samples.
173 if (gpgSamples < MaxGpgSamplesForReset)
174 {
175 reset();
176 return;
177 }
178
179 const GuiderUtils::Vector arc_star = cal.convertToArcseconds(guideStarPosition);
180 const GuiderUtils::Vector arc_reticle = cal.convertToArcseconds(reticlePosition);
181 const GuiderUtils::Vector star_drift = cal.rotateToRaDec(arc_star - arc_reticle);
182
183 double gpgInput = star_drift.x;
184 if (guideStars != nullptr)
185 {
186 double multiStarRADrift, multiStarDECDrift;
187 if (guideStars->getDrift(sqrt(star_drift.x * star_drift.x + star_drift.y * star_drift.y),
188 reticlePosition.x, reticlePosition.y,
189 &multiStarRADrift, &multiStarDECDrift))
190 gpgInput = multiStarRADrift;
191 }
192
193 // Temporarily set the control & prediction gains to 0,
194 // since we won't be guiding while suspended.
195 double predictionGain = gpg->GetPredictionGain();
196 gpg->SetPredictionGain(0.0);
197 double controlGain = gpg->GetControlGain();
198 gpg->SetControlGain(0.0);
199
200 QElapsedTimer gpgTimer;
201 gpgTimer.restart();
202 const double gpgResult = gpg->result(gpgInput, getSNR(guideStars, gpgInput), Options::guideExposure());
203 // Store the updated period length.
204 std::vector<double> gpgParams = gpg->GetGPHyperparameters();
205 Options::setGPGPeriod(gpgParams[PKPeriodLength]);
206 // Not incrementing gpgSamples here, want real samples for that.
207 const double gpgTime = gpgTimer.elapsed();
208 qCDebug(KSTARS_EKOS_GUIDE) << QString("GPG(suspended): elapsed %1s. RA in %2, result: %3")
209 .arg(gpgTime / 1000.0)
210 .arg(gpgInput)
211 .arg(gpgResult);
212 // Reset the gains that were zero'd earlier.
213 gpg->SetPredictionGain(predictionGain);
214 gpg->SetControlGain(controlGain);
215}
216
217bool GPG::computePulse(double raArcsecError, GuideStars *guideStars,
218 int *pulseLength, GuideDirection *pulseDir,
219 const Calibration &cal, Seconds timeStep)
220{
221 if (!Options::gPGEnabled())
222 return false;
223
224 if (fabs(raArcsecError) > MAX_ARCSEC_ERROR)
225 {
226 if (++gpgSkippedSamples > MAX_SKIPPED_SAMPLES)
227 {
228 qCDebug(KSTARS_EKOS_GUIDE) << "Resetting GPG because RA error = "
229 << raArcsecError;
230 reset();
231 }
232 else
233 qCDebug(KSTARS_EKOS_GUIDE) << "Skipping GPG because RA error = "
234 << raArcsecError;
235 return false;
236 }
237 gpgSkippedSamples = 0;
238
239 // I want to start off the gpg on the right foot.
240 // If the initial guiding is messed up by an early focus,
241 // wait until RA has settled down.
242 if (gpgSamples == 0 && fabs(raArcsecError) > 1)
243 {
244 qCDebug(KSTARS_EKOS_GUIDE) << "Delaying GPG startup. RA error = "
245 << raArcsecError;
246 reset();
247 return false;
248 }
249
250 // GPG uses proportional gain and min-move from standard controls. Make sure they're using up-to-date values.
251 gpg->SetControlGain(Options::rAProportionalGain());
252 gpg->SetMinMove(Options::rAMinimumPulseArcSec());
253
254 // GPG input is in RA arcseconds.
255 QElapsedTimer gpgTimer;
256 gpgTimer.restart();
257
258 // Cast back to a raw double
259 auto const rawTime = timeStep.count();
260
261 const double gpgResult = gpg->result(raArcsecError, getSNR(guideStars, raArcsecError), rawTime);
262 const double gpgTime = gpgTimer.elapsed();
263 gpgSamples++;
264
265 // GPG output is in RA arcseconds. pulseLength and pulseDir are set by convertCorrectionToPulseMilliseconds.
266 const double gpgPulse = convertCorrectionToPulseMilliseconds(cal, pulseLength, pulseDir, gpgResult);
267
268 qCDebug(KSTARS_EKOS_GUIDE)
269 << QString("GPG: elapsed %1s. RA in %2, result: %3 * %4 --> %5 : %6ms %7")
270 .arg(gpgTime / 1000.0)
271 .arg(raArcsecError, 0, 'f', 2)
272 .arg(gpgResult, 0, 'f', 2)
273 .arg(cal.raPulseMillisecondsPerArcsecond(), 0, 'f', 1)
274 .arg(gpgPulse, 0, 'f', 1)
275 .arg(*pulseLength)
276 .arg(directionStr(*pulseDir));
277 if (Options::gPGEstimatePeriod())
278 {
279 double period_length = gpg->GetGPHyperparameters()[PKPeriodLength];
280 Options::setGPGPeriod(period_length);
281 }
282 return true;
283}
284
285double GPG::convertCorrectionToPulseMilliseconds(const Calibration &cal, int *pulseLength,
286 GuideDirection *pulseDir, const double gpgResult)
287{
288 const double gpgPulse = gpgResult * cal.raPulseMillisecondsPerArcsecond();
289
290 *pulseDir = gpgPulse > 0 ? RA_DEC_DIR : RA_INC_DIR;
291 *pulseLength = fabs(gpgPulse);
292 return gpgPulse;
293}
294
295bool GPG::darkGuiding(int *pulseLength, GuideDirection *pulseDir, const Calibration &cal,
296 Seconds timeStep)
297{
298 if (!Options::gPGEnabled() || !Options::gPGDarkGuiding())
299 {
300 qCDebug(KSTARS_EKOS_GUIDE) << "dark guiding isn't enabled!";
301 return false;
302 }
303
304 // GPG uses proportional gain and min-move from standard controls. Make sure they're using up-to-date values.
305 gpg->SetControlGain(Options::rAProportionalGain());
306 gpg->SetMinMove(Options::rAMinimumPulseArcSec());
307
308 QElapsedTimer gpgTimer;
309 gpgTimer.restart();
310 const double gpgResult = gpg->deduceResult(timeStep.count());
311 const double gpgTime = gpgTimer.elapsed();
312
313 // GPG output is in RA arcseconds.
314 const double gpgPulse = convertCorrectionToPulseMilliseconds(cal, pulseLength, pulseDir, gpgResult);
315
316 qCDebug(KSTARS_EKOS_GUIDE)
317 << QString("GPG dark guiding: elapsed %1s. RA result: %2 * %3 --> %4 : %5ms %6")
318 .arg(gpgTime / 1000.0)
319 .arg(gpgResult)
320 .arg(cal.raPulseMillisecondsPerArcsecond())
321 .arg(gpgPulse)
322 .arg(*pulseLength)
323 .arg(directionStr(*pulseDir));
324 return true;
325}
This class provides a guiding algorithm for the right ascension axis that learns and predicts the per...
qint64 elapsed() const const
qint64 restart()
QString arg(Args &&... args) const const
Holds all data that is needed for the GP guiding.
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Wed Nov 6 2024 12:13:14 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.