Line data Source code
1 : /** \file clAOLinearPredictor.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Provides a class to manage closed loop gain linear predictor determination.
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : #ifndef clAOLinearPredictor_hpp
9 : #define clAOLinearPredictor_hpp
10 :
11 : #include <cmath>
12 : #include <cstddef>
13 : #include <limits>
14 : #include <vector>
15 :
16 : #include "../../mxlib.hpp"
17 :
18 : #include "../../math/floatUtils.hpp"
19 : #include "../../math/geo.hpp"
20 :
21 : #include "../../sigproc/psdUtils.hpp"
22 :
23 : #include "../../sigproc/autocorrelation.hpp"
24 : #include "../../sigproc/linearPredictor.hpp"
25 :
26 : #include "clGainOpt.hpp"
27 :
28 : namespace mx
29 : {
30 : namespace AO
31 : {
32 : namespace analysis
33 : {
34 :
35 : #define CLAOLP_BREADCRUMB
36 :
37 : // #define CLAOLP_BREADCRUMB std::cerr << __FILE__ << ' ' << __LINE__ << '\n';
38 :
39 : /// Class to manage the calculation of linear predictor coefficients for a closed-loop AO system.
40 : /**
41 : * \tparam _realT the real floating point type in which to do all arithmetic.
42 : *
43 : * \ingroup mxAOAnalytic
44 : */
45 : template <typename _realT>
46 : struct clAOLinearPredictor
47 : {
48 : typedef _realT realT; ///< Floating-point type used for predictor calculations.
49 :
50 : public:
51 : /// Result from one evaluated regularization scale.
52 : struct regResult
53 : {
54 : realT sc; ///< Regularization scale in dB.
55 : realT gopt; ///< Optimum gain at this scale.
56 : realT gmax; ///< Maximum stable gain at this scale.
57 : realT var; ///< Closed-loop variance at the optimum gain.
58 : };
59 :
60 : /// Termination state of the most recent regularization search.
61 : enum class regularizationStatus
62 : {
63 : notRun, ///< No regularization search has been attempted.
64 : converged, ///< The requested precision was reached.
65 : boundaryLimited, ///< The optimum remained on the expanded search boundary.
66 : invalidControls, ///< The configured search controls were invalid.
67 : iterationLimit, ///< The search exhausted its iteration limit.
68 : calculationFailure, ///< Coefficient or gain calculation failed.
69 : };
70 :
71 : /// Diagnostic summary of the most recent regularization search.
72 : struct regularizationReport
73 : {
74 : regularizationStatus status{ regularizationStatus::notRun }; ///< Search termination state.
75 : int iterations{ 0 }; ///< Refinement iterations attempted.
76 : std::size_t evaluations{ 0 }; ///< Regularization scales evaluated.
77 : };
78 :
79 : std::vector<realT> m_PSDtn; ///< Working memory for the regularized PSD
80 :
81 : std::vector<realT> m_psd2s; ///< Working memory for the 2-sided regularized PSD
82 :
83 : std::vector<realT> m_ac; ///< Working memory to hold the autocorrelation.
84 :
85 : sigproc::autocorrelationFromPSD<realT> m_acpsd; ///< Converts the working PSD to an autocorrelation.
86 :
87 : sigproc::linearPredictor<realT> m_lp; ///< Linear predictor used to calculate coefficients.
88 :
89 : realT m_min_var0{ 0 }; ///< Initial minimum variance, with zero requesting initialization.
90 : realT m_min_sc0{ 10 }; ///< Initial minimum regularization scale in dB.
91 : realT m_precision0{ 2 }; ///< Initial regularization scale spacing in dB.
92 : realT m_max_sc0{ 100 }; ///< Initial maximum regularization scale in dB.
93 : realT m_dPrecision{ 3 }; ///< Divisor applied to the spacing during refinement.
94 :
95 : realT m_gmax_lp{ 5 }; ///< The maximum allowable gain for LP.
96 :
97 : // Stopping conditions:
98 : realT m_minPrecision{ 0.001 }; ///< Minimum requested regularization spacing in dB.
99 : int m_maxIts{ 100 }; ///< Maximum number of search refinement iterations.
100 :
101 : int m_extrap{ 1 }; ///< The LP extrapolation length in loop steps. Normally it is 1 step.
102 :
103 : std::vector<regResult> m_regResults; ///< Per-scale telemetry collected when requested.
104 :
105 : regularizationReport m_regularizationReport; ///< Diagnostic summary of the latest search.
106 :
107 : public:
108 : /// Construct a closed-loop linear-predictor calculator with default search controls.
109 20 : clAOLinearPredictor() = default;
110 :
111 : /// Calculate the LP coefficients for a turbulence PSD and a noise PSD.
112 : /** This combines the two PSDs, augments to two-sided, and calls the linearPredictor.calcCoefficients method.
113 : *
114 : * A regularization constant can be added to the PSD as well.
115 : *
116 : * \returns `error_t::noerror` on success, otherwise `error_t::liberr`.
117 : */
118 121 : mx::error_t calcCoefficients( std::vector<realT> &PSDt, /**< [in] the turbulence PSD */
119 : std::vector<realT> &PSDn, /**< [in] the WFS noise PSD */
120 : realT PSDreg, /**< [in] the regularizing constant. Set to 0 to not use. */
121 : int Nc, /**< [in] the number of LP coefficients */
122 : realT condition = 0 /**< [in] the condition number for the SVD. If 0 then
123 : levinson recursion is used. */
124 : )
125 : {
126 : CLAOLP_BREADCRUMB;
127 121 : m_PSDtn.resize( PSDt.size() );
128 :
129 : CLAOLP_BREADCRUMB;
130 7865 : for( size_t i = 0; i < PSDt.size(); ++i )
131 : {
132 7744 : m_PSDtn[i] = PSDt[i] + PSDn[i] + PSDreg;
133 : }
134 :
135 : CLAOLP_BREADCRUMB;
136 121 : sigproc::augment1SidedPSD( m_psd2s, m_PSDtn, 1 );
137 :
138 : CLAOLP_BREADCRUMB;
139 121 : m_ac.resize( m_psd2s.size() );
140 :
141 : CLAOLP_BREADCRUMB;
142 121 : m_acpsd( m_ac, m_psd2s );
143 : CLAOLP_BREADCRUMB;
144 :
145 121 : if( m_lp.calcCoefficients( m_ac, Nc, m_extrap, condition ) != 0 )
146 : {
147 0 : return internal::mxlib_error_report( error_t::liberr, "linearPredictor::calcCoefficients failed" );
148 : }
149 :
150 121 : return error_t::noerror;
151 : }
152 :
153 : /// Worker function for regularizing the PSD for coefficient calculation.
154 : /**
155 : * \tparam telem if true then the results are collected in m_regResults.
156 : *
157 : * On first call (min_var = 0):
158 : * loop over scale factors from min_sc to max_sc (<=) in steps of precision.
159 : *
160 : * On subsequent calls, when min_var and min_sc are passed back in
161 : * loop over scale factors from min_sc-precision to max_sc in steps of
162 : *
163 : * \returns `error_t::noerror` on success, otherwise the coefficient-calculation error.
164 : */
165 : template <bool telem>
166 : mx::error_t
167 6 : _regularizeCoefficients( realT &min_var, /**< [in,out] the minimum variance found; set to 0 on initial call */
168 : realT &min_sc, /**< [in,out] the scale factor at the minimum variance */
169 : realT precision, /**< [in] the step size for the scale factor */
170 : realT max_sc, /**< [in] the maximum scale factor to test */
171 : clGainOpt<realT> &go_lp, /**< [in] the gain optimization object */
172 : std::vector<realT> &PSDt, /**< [in] the turbulence PSD */
173 : std::vector<realT> &PSDn, /**< [in] the WFS noise PSD */
174 : int Nc /**< [in] the number of coefficients */
175 : )
176 : {
177 : CLAOLP_BREADCRUMB;
178 :
179 : realT gmax_lp;
180 : realT gopt_lp;
181 : realT var_lp;
182 :
183 : realT sc0;
184 :
185 6 : if( min_var == 0 ) // first call
186 : {
187 3 : sc0 = min_sc;
188 3 : min_var = std::numeric_limits<realT>::max();
189 : }
190 : else
191 : {
192 3 : sc0 = min_sc - precision * m_dPrecision;
193 : }
194 :
195 : CLAOLP_BREADCRUMB;
196 :
197 : // auto it = std::max_element(std::begin(PSDt), std::end(PSDt));
198 6 : realT psdReg = PSDt[0]; //*it/10;
199 :
200 : CLAOLP_BREADCRUMB;
201 : // Test from sc0 to max_sc in steps of precision
202 : // for( realT sc = sc0; sc <= max_sc; sc += precision )
203 126 : for( realT sc = max_sc; sc >= sc0; sc -= precision )
204 : {
205 : CLAOLP_BREADCRUMB;
206 120 : ++m_regularizationReport.evaluations;
207 120 : error_t rv = calcCoefficients( PSDt, PSDn, psdReg * pow( 10, -sc / 10 ), Nc );
208 120 : if( rv != error_t::noerror )
209 : {
210 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
211 0 : return rv;
212 : }
213 :
214 : CLAOLP_BREADCRUMB;
215 :
216 : CLAOLP_BREADCRUMB;
217 120 : go_lp.a( m_lp.m_c );
218 120 : go_lp.b( m_lp.m_c );
219 :
220 : CLAOLP_BREADCRUMB;
221 120 : rv = go_lp.maxStableGain( gmax_lp );
222 120 : if( rv != error_t::noerror )
223 : {
224 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
225 0 : return rv;
226 : }
227 120 : if( gmax_lp > m_gmax_lp )
228 : {
229 0 : gmax_lp = m_gmax_lp;
230 : }
231 :
232 : CLAOLP_BREADCRUMB;
233 120 : rv = go_lp.optGainOpenLoop( gopt_lp, var_lp, PSDt, PSDn, gmax_lp, false );
234 120 : if( rv != error_t::noerror )
235 : {
236 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
237 0 : return rv;
238 : }
239 :
240 : if( telem )
241 : {
242 120 : m_regResults.push_back( { sc, gopt_lp, gmax_lp, var_lp } );
243 : }
244 :
245 : CLAOLP_BREADCRUMB;
246 120 : if( var_lp < min_var )
247 : {
248 6 : min_var = var_lp;
249 6 : min_sc = sc;
250 : }
251 :
252 : // A jump by a factor of 10 indicates the wall
253 120 : if( var_lp > 10 * min_var )
254 : {
255 0 : return error_t::noerror;
256 : }
257 :
258 : CLAOLP_BREADCRUMB;
259 : }
260 :
261 : CLAOLP_BREADCRUMB;
262 6 : return error_t::noerror;
263 : }
264 :
265 : /// Regularize the PSD and calculate the associated LP coefficients.
266 : /** The PSD is regularized by adding a constant to it. This constant is found by minimizing the variance of the
267 : * residual PSD.
268 : *
269 : * \tparam telem if true then the results are collected in m_regResults
270 : *
271 : * \returns `error_t::noerror` for a converged or boundary-limited search, `error_t::invalidconfig` for invalid
272 : * controls, `error_t::timeout` on iteration exhaustion, or the coefficient-calculation error.
273 : */
274 : template <bool telem = false>
275 : mx::error_t
276 17 : regularizeCoefficients( realT &gmax_lp, /**< [out] the maximum gain calculated for the regularized PSD */
277 : realT &gopt_lp, /**< [out] the optimum gain calculated for the regularized PSD */
278 : realT &var_lp, /**< [out] the variance at the optimum gain */
279 : realT &min_sc, /**< [out] the optimum regularization scale factor */
280 : clGainOpt<realT> &go_lp, /**< [in] the gain optimization object */
281 : std::vector<realT> &PSDt, /**< [in] the turbulence PSD */
282 : std::vector<realT> &PSDn, /**< [in] the WFS noise PSD */
283 : int Nc /**< [in] the number of coefficients */
284 : )
285 : {
286 : CLAOLP_BREADCRUMB;
287 :
288 17 : m_regularizationReport = {};
289 :
290 17 : const realT intervalWidth = m_max_sc0 - m_min_sc0;
291 33 : if( !math::isFinite( m_min_sc0 ) || !math::isFinite( m_max_sc0 ) || !math::isFinite( m_precision0 ) ||
292 13 : !math::isFinite( m_minPrecision ) || !math::isFinite( m_dPrecision ) || m_minPrecision <= 0 ||
293 37 : m_precision0 <= m_minPrecision || intervalWidth <= 0 || m_precision0 > intervalWidth || m_dPrecision <= 1 ||
294 4 : m_maxIts <= 0 )
295 : {
296 14 : m_regularizationReport.status = regularizationStatus::invalidControls;
297 28 : return internal::mxlib_error_report( error_t::invalidconfig,
298 42 : "invalid linear-predictor regularization search controls" );
299 : }
300 :
301 3 : realT min_var = m_min_var0;
302 3 : min_sc = m_min_sc0;
303 3 : realT precision = m_precision0;
304 3 : realT max_sc = m_max_sc0;
305 :
306 : if( telem )
307 : {
308 3 : m_regResults.reserve( m_maxIts * 50 );
309 : }
310 :
311 : CLAOLP_BREADCRUMB;
312 3 : int its = 0;
313 9 : while( precision > m_minPrecision && its < m_maxIts )
314 : {
315 : CLAOLP_BREADCRUMB;
316 6 : const bool firstIteration = its == 0;
317 6 : error_t rv = _regularizeCoefficients<telem>( min_var, min_sc, precision, max_sc, go_lp, PSDt, PSDn, Nc );
318 6 : ++its;
319 6 : m_regularizationReport.iterations = its;
320 6 : if( rv != error_t::noerror )
321 : {
322 0 : return rv;
323 : }
324 :
325 : CLAOLP_BREADCRUMB;
326 6 : if( min_sc == max_sc )
327 : {
328 3 : if( firstIteration )
329 : {
330 3 : min_sc -= precision;
331 3 : max_sc = 200;
332 : }
333 : else
334 : {
335 0 : m_regularizationReport.status = regularizationStatus::boundaryLimited;
336 0 : break;
337 : }
338 : }
339 : else
340 : {
341 3 : max_sc = min_sc + precision;
342 3 : precision /= m_dPrecision;
343 : }
344 : }
345 :
346 3 : if( precision > m_minPrecision && its >= m_maxIts &&
347 2 : m_regularizationReport.status != regularizationStatus::boundaryLimited )
348 : {
349 2 : m_regularizationReport.status = regularizationStatus::iterationLimit;
350 4 : return internal::mxlib_error_report( error_t::timeout,
351 6 : "linear-predictor regularization reached its iteration limit" );
352 : }
353 :
354 1 : if( m_regularizationReport.status != regularizationStatus::boundaryLimited )
355 : {
356 1 : m_regularizationReport.status = regularizationStatus::converged;
357 : }
358 :
359 : CLAOLP_BREADCRUMB;
360 : // Now record final values
361 1 : error_t rv = calcCoefficients( PSDt, PSDn, PSDt[0] * pow( 10, -min_sc / 10 ), Nc );
362 1 : if( rv != error_t::noerror )
363 : {
364 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
365 0 : return rv;
366 : }
367 :
368 : CLAOLP_BREADCRUMB;
369 1 : go_lp.a( m_lp.m_c );
370 1 : go_lp.b( m_lp.m_c );
371 :
372 : CLAOLP_BREADCRUMB;
373 1 : rv = go_lp.maxStableGain( gmax_lp );
374 1 : if( rv != error_t::noerror )
375 : {
376 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
377 0 : return rv;
378 : }
379 :
380 1 : rv = go_lp.optGainOpenLoop( gopt_lp, var_lp, PSDt, PSDn, gmax_lp, false );
381 1 : if( rv != error_t::noerror )
382 : {
383 0 : m_regularizationReport.status = regularizationStatus::calculationFailure;
384 0 : return rv;
385 : }
386 :
387 : CLAOLP_BREADCRUMB;
388 1 : return error_t::noerror;
389 : }
390 :
391 : /// Regularize the PSD and calculate the associated LP coefficients.
392 : /** The PSD is regularized by adding a constant to it. This constant is found by minimizing the variance of the
393 : * residual PSD.
394 : *
395 : * \tparam printout if true then per-scale results are collected in m_regResults.
396 : *
397 : * \returns `error_t::noerror` on success, otherwise the regularization error.
398 : */
399 : template <bool printout = false>
400 : mx::error_t optimizeNc( realT &gmax_lp, /**< [out] maximum gain for the selected predictor */
401 : realT &gopt_lp, /**< [out] optimum gain for the selected predictor */
402 : int &Nc, /**< [out] selected number of coefficients */
403 : realT &var_lp, /**< [out] variance at the optimum gain */
404 : clGainOpt<realT> &go_lp, /**< [in] the gain optimization object */
405 : std::vector<realT> &PSDt, /**< [in] the turbulence PSD */
406 : std::vector<realT> &PSDn, /**< [in] the WFS noise PSD */
407 : int minNc, /**< [in] minimum number of coefficients */
408 : int maxNc /**< [in] maximum number of coefficients */ )
409 : {
410 : realT minVar = std::numeric_limits<realT>::max();
411 :
412 : for( int n = minNc; n <= maxNc; ++n )
413 : {
414 : realT _gmax_lp;
415 : realT _gopt_lp;
416 : realT _var_lp;
417 : realT min_sc;
418 : error_t rv = regularizeCoefficients<printout>( _gmax_lp, _gopt_lp, _var_lp, min_sc, go_lp, PSDt, PSDn, n );
419 : if( rv != error_t::noerror )
420 : {
421 : return rv;
422 : }
423 :
424 : if( _var_lp < minVar )
425 : {
426 : gmax_lp = _gmax_lp;
427 : gopt_lp = _gopt_lp;
428 : var_lp = _var_lp;
429 : Nc = n;
430 :
431 : minVar = var_lp;
432 : }
433 : }
434 :
435 : return error_t::noerror;
436 : }
437 : };
438 :
439 : } // namespace analysis
440 : } // namespace AO
441 : } // namespace mx
442 :
443 : #endif // clAOLinearPredictor_hpp
|
