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 <edgar.klenske@tuebingen.mpg.de>
11 SPDX-FileCopyrightText: Stephan Wenninger <stephan.wenninger@tuebingen.mpg.de>
12 SPDX-FileCopyrightText: Raffi Enficiaud <raffi.enficiaud@tuebingen.mpg.de>
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"
24#include "math_tools.h"
25#include "ekos_guide_debug.h"
26#include <chrono>
27
28enum 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
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
250 /**
251 * Returns the weight of the prediction on the output control value
252 */
254 {
255 auto const period_length = GetGPHyperparameters()[PKPeriodLength];
256 if (get_number_of_measurements() <= 10)
257 {
258 return 0.0;
259 }
260
261 auto current_time = std::chrono::system_clock::now();
262 double delta_measurement_time = std::chrono::duration<double>(current_time - last_time_).count();
263
264 if (parameters.min_periods_for_inference_ * period_length == 0 || delta_measurement_time == 0)
265 {
266 return 0;
267 }
268
269 auto time = std::chrono::duration<double>(current_time - start_time_).count()
270 - (delta_measurement_time / 2.0) // use the midpoint as time stamp
271 + dither_offset_; // correct for the gear time offset from dithering
272 return time / (parameters.min_periods_for_inference_ * period_length);
273 }
274
275
276 GaussianProcessGuider(guide_parameters parameters);
278
279 /**
280 * Calculates the control value based on the current input. 1. The input is
281 * stored, 2. the GP is updated with the new data point, 3. the prediction
282 * is calculated to compensate the gear error and 4. the controller is
283 * calculated, consisting of feedback and prediction parts.
284 */
285 double result(double input, double SNR, double time_step, double prediction_point = -1.0);
286
287 /**
288 * This method provides predictive control if no measurement could be made.
289 * A zero measurement is stored with high uncertainty, and then the GP
290 * prediction is used for control.
291 */
292 double deduceResult(double time_step, double prediction_point = -1.0);
293
294 /**
295 * This method tells the guider that a dither command was issued. The guider
296 * will stop collecting measurements and uses predictions instead, to keep
297 * the FFT and the GP working.
298 */
299 void GuidingDithered(double amt, double rate);
300
301 /**
302 * This method tells the guider that a direct move was
303 * issued. This has the opposite effect of a dither on the dither
304 * offset.
305 */
306 void DirectMoveApplied(double amt, double rate);
307
308 /**
309 * This method tells the guider that dithering is finished. The guider
310 * will resume normal operation.
311 */
312 void GuidingDitherSettleDone(bool success);
313
314 /**
315 * Clears the data from the circular buffer and clears the GP data.
316 */
317 void reset();
318
319 /**
320 * Runs the inference machinery on the GP. Gets the measurement data from
321 * the circular buffer and stores it in Eigen::Vectors. Detrends the data
322 * with linear regression. Calculates the main frequency with an FFT.
323 * Updates the GP accordingly with new data and parameter.
324 */
325 void UpdateGP(double prediction_point = std::numeric_limits<double>::quiet_NaN());
326
327 /**
328 * Does filtering and sets the period length of the GPGuider.
329 */
331
332 data_point &get_last_point() const
333 {
334 return circular_buffer_data_[circular_buffer_data_.size() - 1];
335 }
336
337 data_point &get_second_last_point() const
338 {
339 return circular_buffer_data_[circular_buffer_data_.size() - 2];
340 }
341
342 size_t get_number_of_measurements() const
343 {
344 return circular_buffer_data_.size();
345 }
346
347 void add_one_point()
348 {
349 circular_buffer_data_.push_front(data_point());
350 }
351
352 /**
353 * This method is needed for automated testing. It can inject data points.
354 */
355 void inject_data_point(double timestamp, double input, double SNR, double control);
356
357 /**
358 * Takes timestamps, measurements and SNRs and returns them regularized in a matrix.
359 */
360 Eigen::MatrixXd regularize_dataset(const Eigen::VectorXd &timestamps, const Eigen::VectorXd &gear_error,
361 const Eigen::VectorXd &variances);
362
363 /**
364 * Saves the GP data to a csv file for external analysis. Expensive!
365 */
366 void save_gp_data() const;
367
368 /**
369 * Sets the learning rate. Useful for disabling it for testing.
370 */
371 void SetLearningRate(double learning_rate);
372};
373
374//
375// GPDebug abstract interface to allow logging to the PHD2 Debug Log when
376// GaussianProcessGuider is built in the context of the PHD2 application
377//
378// Add code like this to record debug info in the PHD2 debug log (with newline appended)
379//
380// GPDebug->Log("input: %.2f SNR: %.1f time_step: %.1f", input, SNR, time_step);
381//
382// Outside of PHD2, like in the test framework, these calls will not produce any output
383//
384class GPDebug
385{
386 public:
387 static void SetGPDebug(GPDebug *logger);
388 virtual ~GPDebug();
389 virtual void Log(const char *format, ...) = 0;
390};
391
392extern class GPDebug *GPDebug;
393
394#endif // GAUSSIAN_PROCESS_GUIDER
This class provides a guiding algorithm for the right ascension axis that learns and predicts the per...
void UpdateGP(double prediction_point=std::numeric_limits< double >::quiet_NaN())
Runs the inference machinery on the GP.
void SetLearningRate(double learning_rate)
Sets the learning rate.
void GuidingDitherSettleDone(bool success)
This method tells the guider that dithering is finished.
double getPredictionContribution()
Returns the weight of the prediction on the output control value.
void reset()
Clears the data from the circular buffer and clears the GP data.
void save_gp_data() const
Saves the GP data to a csv file for external analysis.
void GuidingDithered(double amt, double rate)
This method tells the guider that a dither command was issued.
void DirectMoveApplied(double amt, double rate)
This method tells the guider that a direct move was issued.
void UpdatePeriodLength(double period_length)
Does filtering and sets the period length of the GPGuider.
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.
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.
double result(double input, double SNR, double time_step, double prediction_point=-1.0)
Calculates the control value based on the current input.
Implementation of Open Astronomy Log (OAL) XML specifications to record observation logs.
The file holds the covariance functions that can be used with the GP class.
The GP class implements the Gaussian Process functionality.
Provides mathematical tools needed for the Gaussian process toolbox.
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 Fri May 17 2024 11:48:26 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.