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 namespace
17 {
18 
19 // Don't have GPG model large errors. Instantaneous large errors are probably not due to
20 // periodic error, but rather one-offs.
21 constexpr double MAX_ARCSEC_ERROR = 5.0;
22 
23 // If we skip this many samples, reset gpg.
24 constexpr int MAX_SKIPPED_SAMPLES = 4;
25 
26 // Fills parameters from the KStars options.
27 void getGPGParameters(GaussianProcessGuider::guide_parameters *parameters)
28 {
29  // Parameters from standard options.
30  parameters->control_gain_ = Options::rAProportionalGain();
31  // We need a calibration to determine min move. This is re-set in computePulse() below.
32  parameters->min_move_ = 0.25;
33 
34  // Parameters from GPG-specific options.
35  parameters->min_periods_for_inference_ = Options::gPGMinPeriodsForInference();
36  parameters->SE0KLengthScale_ = Options::gPGSE0KLengthScale();
37  parameters->SE0KSignalVariance_ = Options::gPGSE0KSignalVariance();
38  parameters->PKLengthScale_ = Options::gPGPKLengthScale();
39  parameters->PKPeriodLength_ = Options::gPGPeriod();
40  parameters->PKSignalVariance_ = Options::gPGPKSignalVariance();
41  parameters->SE1KLengthScale_ = Options::gPGSE1KLengthScale();
42  parameters->SE1KSignalVariance_ = Options::gPGSE1KSignalVariance();
43  parameters->min_periods_for_period_estimation_ = Options::gPGMinPeriodsForPeriodEstimate();
44  parameters->points_for_approximation_ = Options::gPGPointsForApproximation();
45  parameters->prediction_gain_ = Options::gPGpWeight();
46  parameters->compute_period_ = Options::gPGEstimatePeriod();
47 }
48 
49 // Returns the SNR returned by guideStars, or if guideStars is null (e.g. we aren't
50 // running SEP MultiStar, then just 50, as this isn't a critical parameter.
51 double getSNR(const GuideStars *guideStars, const double raArcsecError)
52 {
53  if (fabs(raArcsecError) > MAX_ARCSEC_ERROR)
54  return 1.0;
55 
56  // If no SEP MultiStar, just fudge the SNR. Won't be a big deal.
57  if (guideStars == nullptr)
58  return 50.0;
59  else
60  return guideStars->getGuideStarSNR();
61 }
62 
63 } // namespace
64 
65 void GPG::updateParameters()
66 {
67  // Parameters would be set when the gpg stars up.
68  if (gpg.get() == nullptr) return;
69 
71  getGPGParameters(&parameters);
72  gpg->SetControlGain(parameters.control_gain_);
73  gpg->SetPeriodLengthsInference(parameters.min_periods_for_inference_);
74  gpg->SetMinMove(parameters.min_move_);
75  gpg->SetPeriodLengthsPeriodEstimation(parameters.min_periods_for_period_estimation_);
76  gpg->SetNumPointsForApproximation(parameters.points_for_approximation_);
77  gpg->SetPredictionGain(parameters.prediction_gain_);
78  gpg->SetBoolComputePeriod(parameters.compute_period_);
79 
80  // The GPG header really should be in a namespace so NumParameters
81  // is not in a global namespace.
82  std::vector<double> hyperparameters(NumParameters);
83  hyperparameters[SE0KLengthScale] = parameters.SE0KLengthScale_;
84  hyperparameters[SE0KSignalVariance] = parameters.SE0KSignalVariance_;
85  hyperparameters[PKLengthScale] = parameters.PKLengthScale_;
86  hyperparameters[PKSignalVariance] = parameters.PKSignalVariance_;
87  hyperparameters[SE1KLengthScale] = parameters.SE1KLengthScale_;
88  hyperparameters[SE1KSignalVariance] = parameters.SE1KSignalVariance_;
89  hyperparameters[PKPeriodLength] = parameters.PKPeriodLength_;
90  gpg->SetGPHyperparameters(hyperparameters);
91 }
92 
93 
94 GPG::GPG()
95 {
97  getGPGParameters(&parameters);
98  gpg.reset(new GaussianProcessGuider(parameters));
99  reset();
100 }
101 
102 void GPG::reset()
103 {
104  gpgSamples = 0;
105  gpgSkippedSamples = 0;
106  gpg->reset();
107  qCDebug(KSTARS_EKOS_GUIDE) << "Resetting GPG";
108 }
109 
110 void GPG::startDithering(double dx, double dy, const Calibration &cal)
111 {
112  // convert the x and y offsets to RA and DEC offsets.
113  double raPixels, decPixels;
114  cal.rotateToRaDec(dx, dy, &raPixels, &decPixels);
115 
116  // The 1000 in the denominator below converts gear-milliseconds into gear-seconds,
117  // which is the unit GPG wants. Since calibration uses arcseconds, we convert using
118  // arcseconds/pixel, though that's inexact when the pixels aren't square.
119  // amount = (pixels * ms/pixel) / 1000;
120  const double amount = raPixels * cal.raPulseMillisecondsPerArcsecond() * cal.xArcsecondsPerPixel() / 1000.0;
121 
122  qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither started. Gear-seconds =" << amount << "x,y was: " << dx << dy;
123  gpg->GuidingDithered(amount, 1.0);
124 }
125 
126 void GPG::ditheringSettled(bool success)
127 {
128  gpg->GuidingDitherSettleDone(success);
129  if (success)
130  qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither done (success)";
131  else
132  qCDebug(KSTARS_EKOS_GUIDE) << "GPG Dither done (failed)";
133 }
134 
135 void GPG::suspended(const GuiderUtils::Vector &guideStarPosition,
136  const GuiderUtils::Vector &reticlePosition,
137  GuideStars *guideStars,
138  const Calibration &cal)
139 {
140  constexpr int MaxGpgSamplesForReset = 25;
141  // We just reset the gpg if there's not enough samples to make
142  // it worth our while to keep it going. We're probably worse off
143  // to start it and then feed it bad samples.
144  if (gpgSamples < MaxGpgSamplesForReset)
145  {
146  reset();
147  return;
148  }
149 
150  const GuiderUtils::Vector arc_star = cal.convertToArcseconds(guideStarPosition);
151  const GuiderUtils::Vector arc_reticle = cal.convertToArcseconds(reticlePosition);
152  const GuiderUtils::Vector star_drift = cal.rotateToRaDec(arc_star - arc_reticle);
153 
154  double gpgInput = star_drift.x;
155  if (guideStars != nullptr)
156  {
157  double multiStarRADrift, multiStarDECDrift;
158  if (guideStars->getDrift(sqrt(star_drift.x * star_drift.x + star_drift.y * star_drift.y),
159  reticlePosition.x, reticlePosition.y,
160  &multiStarRADrift, &multiStarDECDrift))
161  gpgInput = multiStarRADrift;
162  }
163 
164  // Temporarily set the control & prediction gains to 0,
165  // since we won't be guiding while suspended.
166  double predictionGain = gpg->GetPredictionGain();
167  gpg->SetPredictionGain(0.0);
168  double controlGain = gpg->GetControlGain();
169  gpg->SetControlGain(0.0);
170 
171  QElapsedTimer gpgTimer;
172  gpgTimer.restart();
173  const double gpgResult = gpg->result(gpgInput, getSNR(guideStars, gpgInput), Options::guideExposure());
174  // Store the updated period length.
175  std::vector<double> gpgParams = gpg->GetGPHyperparameters();
176  Options::setGPGPeriod(gpgParams[PKPeriodLength]);
177  // Not incrementing gpgSamples here, want real samples for that.
178  const double gpgTime = gpgTimer.elapsed();
179  qCDebug(KSTARS_EKOS_GUIDE) << QString("GPG(suspended): elapsed %1s. RA in %2, result: %3")
180  .arg(gpgTime / 1000.0)
181  .arg(gpgInput)
182  .arg(gpgResult);
183  // Reset the gains that were zero'd earlier.
184  gpg->SetPredictionGain(predictionGain);
185  gpg->SetControlGain(controlGain);
186 }
187 
188 bool GPG::computePulse(double raArcsecError, GuideStars *guideStars,
189  int *pulseLength, GuideDirection *pulseDir,
190  const Calibration &cal)
191 {
192  if (!Options::gPGEnabled())
193  return false;
194 
195  if (fabs(raArcsecError) > MAX_ARCSEC_ERROR)
196  {
197  if (++gpgSkippedSamples > MAX_SKIPPED_SAMPLES)
198  {
199  qCDebug(KSTARS_EKOS_GUIDE) << "Resetting GPG because RA error = "
200  << raArcsecError;
201  reset();
202  }
203  else
204  qCDebug(KSTARS_EKOS_GUIDE) << "Skipping GPG because RA error = "
205  << raArcsecError;
206  return false;
207  }
208  gpgSkippedSamples = 0;
209 
210  // I want to start off the gpg on the right foot.
211  // If the initial guiding is messed up by an early focus,
212  // wait until RA has settled down.
213  if (gpgSamples == 0 && fabs(raArcsecError) > 1)
214  {
215  qCDebug(KSTARS_EKOS_GUIDE) << "Delaying GPG startup. RA error = "
216  << raArcsecError;
217  reset();
218  return false;
219  }
220 
221  // GPG uses proportional gain and min-move from standard controls. Make sure they're using up-to-date values.
222  gpg->SetControlGain(Options::rAProportionalGain());
223  gpg->SetMinMove(Options::rAMinimumPulseArcSec());
224 
225  // GPG input is in RA arcseconds.
226  QElapsedTimer gpgTimer;
227  gpgTimer.restart();
228  const double gpgResult = gpg->result(raArcsecError, getSNR(guideStars, raArcsecError), Options::guideExposure());
229  const double gpgTime = gpgTimer.elapsed();
230  gpgSamples++;
231 
232  // GPG output is in RA arcseconds.
233  const double gpgPulse = gpgResult * cal.raPulseMillisecondsPerArcsecond();
234 
235  *pulseDir = gpgPulse > 0 ? RA_DEC_DIR : RA_INC_DIR;
236  *pulseLength = fabs(gpgPulse);
237 
238  if (*pulseDir == RA_DEC_DIR)
239  qCDebug(KSTARS_EKOS_GUIDE) << "pulse_length [RA] =" << *pulseLength << "Direction : Decrease";
240  else
241  qCDebug(KSTARS_EKOS_GUIDE) << "pulse_length [RA] =" << *pulseLength << "Direction : Increase";
242 
243  qCDebug(KSTARS_EKOS_GUIDE)
244  << QString("GPG: elapsed %1s. RA in %2, result: %3 * %4 --> %5 : len %6 dir %7")
245  .arg(gpgTime / 1000.0)
246  .arg(raArcsecError)
247  .arg(gpgResult)
248  .arg(cal.raPulseMillisecondsPerArcsecond())
249  .arg(gpgPulse)
250  .arg(*pulseLength)
251  .arg(*pulseDir);
252  if (Options::gPGEstimatePeriod())
253  {
254  double period_length = gpg->GetGPHyperparameters()[PKPeriodLength];
255  Options::setGPGPeriod(period_length);
256  }
257  return true;
258 }
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-2022 The KDE developers.
Generated on Fri Aug 19 2022 03:57:51 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.