Kstars

gpg.cpp
1 /*
2  SPDX-FileCopyrightText: 2020 Hy Murveit <[email protected]>
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 
18 namespace
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.
23 constexpr double MAX_ARCSEC_ERROR = 5.0;
24 
25 // If we skip this many samples, reset gpg.
26 constexpr int MAX_SKIPPED_SAMPLES = 4;
27 
28 // Fills parameters from the KStars options.
29 void 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.
53 double 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 
65 const 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 
84 double 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 
94 void 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 
123 GPG::GPG()
124 {
126  getGPGParameters(&parameters);
127  gpg.reset(new GaussianProcessGuider(parameters));
128  reset();
129 }
130 
131 void GPG::reset()
132 {
133  gpgSamples = 0;
134  gpgSkippedSamples = 0;
135  gpg->reset();
136  qCDebug(KSTARS_EKOS_GUIDE) << "Resetting GPG";
137 }
138 
139 void 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 
155 void 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 
164 void 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 
217 bool 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 
285 double 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 
295 bool 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 }
qint64 restart()
This class provides a guiding algorithm for the right ascension axis that learns and predicts the per...
qint64 elapsed() const const
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
KGuiItem reset()
Holds all data that is needed for the GP guiding.
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Tue Sep 26 2023 03:55:46 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.