Kstars

gaussian_process_guider.h
1 /*
2  PHD2 Guiding
3 
4  @file
5  @date 2014-2017
6  @copyright Max Planck Society
7 
8  @brief Provides a Gaussian process based guiding algorithm
9 
10  SPDX-FileCopyrightText: Edgar D. Klenske <[email protected]>
11  SPDX-FileCopyrightText: Stephan Wenninger <[email protected]>
12  SPDX-FileCopyrightText: Raffi Enficiaud <[email protected]>
13  SPDX-License-Identifier: BSD-3-Clause
14 */
15 
16 #ifndef GAUSSIAN_PROCESS_GUIDER
17 #define GAUSSIAN_PROCESS_GUIDER
18 
19 // Removed the dependence on phd2
20 // #include "circbuf.h"
21 
22 #include "gaussian_process.h"
23 #include "covariance_functions.h"
24 #include "math_tools.h"
25 
26 #include <chrono>
27 
28 enum Hyperparameters
29 {
30  SE0KLengthScale,
31  SE0KSignalVariance,
32  PKLengthScale,
33  PKSignalVariance,
34  SE1KLengthScale,
35  SE1KSignalVariance,
36  PKPeriodLength,
37  NumParameters // this represents the number of elements in the enum
38 };
39 
40 /**
41  * This class provides a guiding algorithm for the right ascension axis that
42  * learns and predicts the periodic gear error with a Gaussian process. This
43  * prediction helps reducing periodic error components in the residual tracking
44  * error. Further it is able to perform tracking without measurement to increase
45  * robustness of the overall guiding system.
46  */
48 {
49  // HY: Added this subclass to remove dependency on phd2.
50  // Other than this class, this is the original GPG code.
51  // This implements a basix circular buffer that can hold upto size items.
52  // Cannot resize, cannot remove items. Can just insert and randomly access
53  // the last "size" items in the order inserted, or clear the entire buffer.
54  template <typename T>
55  class circular_buffer
56  {
57  public:
58  circular_buffer(int size) : size_(size), buff(new T[size]) {}
59  ~circular_buffer()
60  {
61  delete[] buff;
62  }
63 
64  T &operator[](int index) const
65  {
66  assert(index >= 0 && index < size_ && index < num_items);
67  int location = (start + index ) % size_;
68  return buff[location];
69  }
70 
71  int size() const
72  {
73  return num_items;
74  }
75 
76  void clear()
77  {
78  start = 0;
79  end = 0;
80  num_items = 0;
81  }
82 
83  void push_front(T item)
84  {
85  buff[end] = item;
86  num_items++;
87  if (num_items >= size_) num_items = size_;
88  end++;
89  if (end >= size_) end = 0;
90  }
91 
92  private:
93  int size_ = 0; // Max length of the queue.
94  T *buff; // Data storage.
95  int start = 0; // Index of oldest entry in the buffer..
96  int end = 0; // Index of the next place to put data into the buffer.
97  int num_items = 0; // Number of items currently in the circular buffer.
98  };
99 
100  public:
101 
102  struct data_point
103  {
104  double timestamp;
105  double measurement; // current pointing error
106  double variance; // current measurement variance
107  double control; // control action
108  };
109 
110  /**
111  * Holds all data that is needed for the GP guiding.
112  */
114  {
115  double control_gain_;
116  double min_move_;
117  double prediction_gain_;
118 
119  double min_periods_for_inference_;
120  double min_periods_for_period_estimation_;
121 
122  int points_for_approximation_;
123 
124  bool compute_period_;
125 
126  double SE0KLengthScale_;
127  double SE0KSignalVariance_;
128  double PKLengthScale_;
129  double PKSignalVariance_;
130  double SE1KLengthScale_;
131  double SE1KSignalVariance_;
132  double PKPeriodLength_;
133 
134  guide_parameters() :
135  control_gain_(0.0),
136  min_move_(0.0),
137  prediction_gain_(0.0),
138  min_periods_for_inference_(0.0),
139  min_periods_for_period_estimation_(0.0),
140  points_for_approximation_(0),
141  compute_period_(false),
142  SE0KLengthScale_(0.0),
143  SE0KSignalVariance_(0.0),
144  PKLengthScale_(0.0),
145  PKSignalVariance_(0.0),
146  SE1KLengthScale_(0.0),
147  SE1KSignalVariance_(0.0),
148  PKPeriodLength_(0.0)
149  {
150  }
151 
152  };
153 
154  private:
155 
156  std::chrono::system_clock::time_point start_time_; // reference time
157  std::chrono::system_clock::time_point last_time_;
158 
159  double control_signal_;
160  double prediction_;
161  double last_prediction_end_;
162 
163  int dither_steps_;
164  bool dithering_active_;
165 
166  //! the dither offset collects the correction in gear time from dithering
167  double dither_offset_;
168 
169  circular_buffer<data_point> circular_buffer_data_;
170 
171  covariance_functions::PeriodicSquareExponential2 covariance_function_; // for inference
172  covariance_functions::PeriodicSquareExponential output_covariance_function_; // for prediction
173  GP gp_;
174 
175  /**
176  * Learning rate for smooth parameter adaptation.
177  */
178  double learning_rate_;
179 
180  /**
181  * Guiding parameters of this instance.
182  */
183  guide_parameters parameters;
184 
185  /**
186  * Stores the current time and creates a timestamp for the GP.
187  */
188  void SetTimestamp();
189 
190  /**
191  * Stores the measurement, SNR and resets last_prediction_end_.
192  */
193  void HandleGuiding(double input, double SNR);
194 
195  /**
196  * Stores a zero as blind "measurement" with high variance.
197  */
198  void HandleDarkGuiding();
199 
200  /**
201  * Stores the control value.
202  */
203  void HandleControls(double control_input);
204 
205  /**
206  * Calculates the noise from the reported SNR value according to an
207  * empirically justified equation and stores it.
208  */
209  double CalculateVariance(double SNR);
210 
211  /**
212  * Estimates the main period length for a given dataset.
213  */
214  double EstimatePeriodLength(const Eigen::VectorXd &time, const Eigen::VectorXd &data);
215 
216  /**
217  * Calculates the difference in gear error for the time between the last
218  * prediction point and the current prediction point, which lies one
219  * exposure length in the future.
220  */
221  double PredictGearError(double prediction_location);
222 
223 
224 
225  public:
226  double GetControlGain() const;
227  bool SetControlGain(double control_gain);
228 
229  double GetMinMove() const;
230  bool SetMinMove(double min_move);
231 
232  double GetPeriodLengthsInference() const;
233  bool SetPeriodLengthsInference(double num_periods);
234 
235  double GetPeriodLengthsPeriodEstimation() const;
236  bool SetPeriodLengthsPeriodEstimation(double num_periods);
237 
238  int GetNumPointsForApproximation() const;
239  bool SetNumPointsForApproximation(int num_points);
240 
241  bool GetBoolComputePeriod() const;
242  bool SetBoolComputePeriod(bool active);
243 
244  std::vector<double> GetGPHyperparameters() const;
245  bool SetGPHyperparameters(const std::vector<double> &hyperparameters);
246 
247  double GetPredictionGain() const;
248  bool SetPredictionGain(double);
249 
252 
253  /**
254  * Calculates the control value based on the current input. 1. The input is
255  * stored, 2. the GP is updated with the new data point, 3. the prediction
256  * is calculated to compensate the gear error and 4. the controller is
257  * calculated, consisting of feedback and prediction parts.
258  */
259  double result(double input, double SNR, double time_step, double prediction_point = -1.0);
260 
261  /**
262  * This method provides predictive control if no measurement could be made.
263  * A zero measurement is stored with high uncertainty, and then the GP
264  * prediction is used for control.
265  */
266  double deduceResult(double time_step, double prediction_point = -1.0);
267 
268  /**
269  * This method tells the guider that a dither command was issued. The guider
270  * will stop collecting measurements and uses predictions instead, to keep
271  * the FFT and the GP working.
272  */
273  void GuidingDithered(double amt, double rate);
274 
275  /**
276  * This method tells the guider that a direct move was
277  * issued. This has the opposite effect of a dither on the dither
278  * offset.
279  */
280  void DirectMoveApplied(double amt, double rate);
281 
282  /**
283  * This method tells the guider that dithering is finished. The guider
284  * will resume normal operation.
285  */
286  void GuidingDitherSettleDone(bool success);
287 
288  /**
289  * Clears the data from the circular buffer and clears the GP data.
290  */
291  void reset();
292 
293  /**
294  * Runs the inference machinery on the GP. Gets the measurement data from
295  * the circular buffer and stores it in Eigen::Vectors. Detrends the data
296  * with linear regression. Calculates the main frequency with an FFT.
297  * Updates the GP accordingly with new data and parameter.
298  */
299  void UpdateGP(double prediction_point = std::numeric_limits<double>::quiet_NaN());
300 
301  /**
302  * Does filtering and sets the period length of the GPGuider.
303  */
304  void UpdatePeriodLength(double period_length);
305 
306  data_point &get_last_point() const
307  {
308  return circular_buffer_data_[circular_buffer_data_.size() - 1];
309  }
310 
311  data_point &get_second_last_point() const
312  {
313  return circular_buffer_data_[circular_buffer_data_.size() - 2];
314  }
315 
316  size_t get_number_of_measurements() const
317  {
318  return circular_buffer_data_.size();
319  }
320 
321  void add_one_point()
322  {
323  circular_buffer_data_.push_front(data_point());
324  }
325 
326  /**
327  * This method is needed for automated testing. It can inject data points.
328  */
329  void inject_data_point(double timestamp, double input, double SNR, double control);
330 
331  /**
332  * Takes timestamps, measurements and SNRs and returns them regularized in a matrix.
333  */
334  Eigen::MatrixXd regularize_dataset(const Eigen::VectorXd &timestamps, const Eigen::VectorXd &gear_error,
335  const Eigen::VectorXd &variances);
336 
337  /**
338  * Saves the GP data to a csv file for external analysis. Expensive!
339  */
340  void save_gp_data() const;
341 
342  /**
343  * Sets the learning rate. Useful for disabling it for testing.
344  */
345  void SetLearningRate(double learning_rate);
346 };
347 
348 //
349 // GPDebug abstract interface to allow logging to the PHD2 Debug Log when
350 // GaussianProcessGuider is built in the context of the PHD2 application
351 //
352 // Add code like this to record debug info in the PHD2 debug log (with newline appended)
353 //
354 // GPDebug->Log("input: %.2f SNR: %.1f time_step: %.1f", input, SNR, time_step);
355 //
356 // Outside of PHD2, like in the test framework, these calls will not produce any output
357 //
358 class GPDebug
359 {
360  public:
361  static void SetGPDebug(GPDebug *logger);
362  virtual ~GPDebug();
363  virtual void Log(const char *format, ...) = 0;
364 };
365 
366 extern class GPDebug *GPDebug;
367 
368 #endif // GAUSSIAN_PROCESS_GUIDER
The GP class implements the Gaussian Process functionality.
void DirectMoveApplied(double amt, double rate)
This method tells the guider that a direct move was issued.
double result(double input, double SNR, double time_step, double prediction_point=-1.0)
Calculates the control value based on the current input.
The file holds the covariance functions that can be used with the GP class.
void reset()
Clears the data from the circular buffer and clears the GP data.
void UpdatePeriodLength(double period_length)
Does filtering and sets the period length of the GPGuider.
void UpdateGP(double prediction_point=std::numeric_limits< double >::quiet_NaN())
Runs the inference machinery on the GP.
double deduceResult(double time_step, double prediction_point=-1.0)
This method provides predictive control if no measurement could be made.
void inject_data_point(double timestamp, double input, double SNR, double control)
This method is needed for automated testing.
void GuidingDithered(double amt, double rate)
This method tells the guider that a dither command was issued.
This class provides a guiding algorithm for the right ascension axis that learns and predicts the per...
void SetLearningRate(double learning_rate)
Sets the learning rate.
Eigen::MatrixXd regularize_dataset(const Eigen::VectorXd &timestamps, const Eigen::VectorXd &gear_error, const Eigen::VectorXd &variances)
Takes timestamps, measurements and SNRs and returns them regularized in a matrix.
Holds all data that is needed for the GP guiding.
void save_gp_data() const
Saves the GP data to a csv file for external analysis.
void GuidingDitherSettleDone(bool success)
This method tells the guider that dithering is finished.
Provides mathematical tools needed for the Gaussian process toolbox.
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:54 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.