From 2aacb6f175c5d38467a67db0326bc2ec5fced4b0 Mon Sep 17 00:00:00 2001 From: Jared Males Date: Sun, 2 Aug 2026 10:53:51 -0700 Subject: [PATCH 1/2] Harden linear predictor regularization controls --- include/ao/analysis/clAOLinearPredictor.hpp | 237 ++++++++++++------ include/ao/analysis/fourierTemporalPSD.hpp | 22 +- include/improc/imageUtils.hpp | 27 +- include/math/CMakeLists.txt | 1 + include/math/floatUtils.hpp | 95 +++++++ include/math/math.hpp | 1 + tests/CMakeLists.txt | 2 + tests/Makefile | 2 + .../ao/analysis/clAOLinearPredictor_test.cpp | 209 +++++++++++++++ tests/include/improc/imageUtils_test.cpp | 16 +- tests/include/math/floatUtils_test.cpp | 31 +++ 11 files changed, 535 insertions(+), 108 deletions(-) create mode 100644 include/math/floatUtils.hpp create mode 100644 tests/include/ao/analysis/clAOLinearPredictor_test.cpp create mode 100644 tests/include/math/floatUtils_test.cpp diff --git a/include/ao/analysis/clAOLinearPredictor.hpp b/include/ao/analysis/clAOLinearPredictor.hpp index 852bcffba..a2594ae60 100644 --- a/include/ao/analysis/clAOLinearPredictor.hpp +++ b/include/ao/analysis/clAOLinearPredictor.hpp @@ -8,8 +8,14 @@ #ifndef clAOLinearPredictor_hpp #define clAOLinearPredictor_hpp +#include +#include +#include #include +#include "../../mxlib.hpp" + +#include "../../math/floatUtils.hpp" #include "../../math/geo.hpp" #include "../../sigproc/psdUtils.hpp" @@ -28,8 +34,7 @@ namespace analysis #define CLAOLP_BREADCRUMB -//#define CLAOLP_BREADCRUMB std::cerr << __FILE__ << ' ' << __LINE__ << '\n'; - +// #define CLAOLP_BREADCRUMB std::cerr << __FILE__ << ' ' << __LINE__ << '\n'; /// Class to manage the calculation of linear predictor coefficients for a closed-loop AO system. /** @@ -40,61 +45,83 @@ namespace analysis template struct clAOLinearPredictor { - typedef _realT realT; - -public: + typedef _realT realT; ///< Floating-point type used for predictor calculations. + public: + /// Result from one evaluated regularization scale. struct regResult { - realT sc; - realT gopt; - realT gmax; - realT var; + realT sc; ///< Regularization scale in dB. + realT gopt; ///< Optimum gain at this scale. + realT gmax; ///< Maximum stable gain at this scale. + realT var; ///< Closed-loop variance at the optimum gain. + }; + + /// Termination state of the most recent regularization search. + enum class regularizationStatus + { + notRun, ///< No regularization search has been attempted. + converged, ///< The requested precision was reached. + boundaryLimited, ///< The optimum remained on the expanded search boundary. + invalidControls, ///< The configured search controls were invalid. + iterationLimit, ///< The search exhausted its iteration limit. + calculationFailure, ///< Coefficient or gain calculation failed. + }; + + /// Diagnostic summary of the most recent regularization search. + struct regularizationReport + { + regularizationStatus status{ regularizationStatus::notRun }; ///< Search termination state. + int iterations{ 0 }; ///< Refinement iterations attempted. + std::size_t evaluations{ 0 }; ///< Regularization scales evaluated. }; - std::vector m_PSDtn; ///< Working memory for the regularized PSD + std::vector m_PSDtn; ///< Working memory for the regularized PSD - std::vector m_psd2s; ///< Working memory for the 2-sided regularized PSD + std::vector m_psd2s; ///< Working memory for the 2-sided regularized PSD - std::vector m_ac; ///< Working memory to hold the autocorrelation. + std::vector m_ac; ///< Working memory to hold the autocorrelation. - sigproc::autocorrelationFromPSD m_acpsd; + sigproc::autocorrelationFromPSD m_acpsd; ///< Converts the working PSD to an autocorrelation. - sigproc::linearPredictor m_lp; + sigproc::linearPredictor m_lp; ///< Linear predictor used to calculate coefficients. - realT m_min_var0{ 0 }; - realT m_min_sc0{ 10 }; - realT m_precision0{ 2 }; - realT m_max_sc0{ 100 }; - realT m_dPrecision{ 3 }; + realT m_min_var0{ 0 }; ///< Initial minimum variance, with zero requesting initialization. + realT m_min_sc0{ 10 }; ///< Initial minimum regularization scale in dB. + realT m_precision0{ 2 }; ///< Initial regularization scale spacing in dB. + realT m_max_sc0{ 100 }; ///< Initial maximum regularization scale in dB. + realT m_dPrecision{ 3 }; ///< Divisor applied to the spacing during refinement. - realT m_gmax_lp{ 5 }; ///< The maximum allowable gain for LP. + realT m_gmax_lp{ 5 }; ///< The maximum allowable gain for LP. // Stopping conditions: - realT m_minPrecision{ 0.001 }; - int m_maxIts{ 100 }; + realT m_minPrecision{ 0.001 }; ///< Minimum requested regularization spacing in dB. + int m_maxIts{ 100 }; ///< Maximum number of search refinement iterations. + + int m_extrap{ 1 }; ///< The LP extrapolation length in loop steps. Normally it is 1 step. - int m_extrap {1}; ///< The LP extrapolation length in loop steps. Normally it is 1 step. + std::vector m_regResults; ///< Per-scale telemetry collected when requested. - std::vector m_regResults; -public: + regularizationReport m_regularizationReport; ///< Diagnostic summary of the latest search. - clAOLinearPredictor() - {} + public: + /// Construct a closed-loop linear-predictor calculator with default search controls. + clAOLinearPredictor() = default; /// Calculate the LP coefficients for a turbulence PSD and a noise PSD. /** This combines the two PSDs, augments to two-sided, and calls the linearPredictor.calcCoefficients method. * * A regularization constant can be added to the PSD as well. * + * \returns `error_t::noerror` on success, otherwise `error_t::liberr`. */ - int calcCoefficients( std::vector &PSDt, ///< [in] the turbulence PSD - std::vector &PSDn, ///< [in] the WFS noise PSD - realT PSDreg, ///< [in] the regularizing constant. Set to 0 to not use. - int Nc, ///< [in] the number of LP coefficients. - realT condition = 0 /**< [in] the condition number for the SVD. If 0 then - levinson recursion is used. */ - ) + mx::error_t calcCoefficients( std::vector &PSDt, /**< [in] the turbulence PSD */ + std::vector &PSDn, /**< [in] the WFS noise PSD */ + realT PSDreg, /**< [in] the regularizing constant. Set to 0 to not use. */ + int Nc, /**< [in] the number of LP coefficients */ + realT condition = 0 /**< [in] the condition number for the SVD. If 0 then + levinson recursion is used. */ + ) { CLAOLP_BREADCRUMB; m_PSDtn.resize( PSDt.size() ); @@ -115,8 +142,12 @@ struct clAOLinearPredictor m_acpsd( m_ac, m_psd2s ); CLAOLP_BREADCRUMB; - return m_lp.calcCoefficients( m_ac, Nc, m_extrap , condition ); + if( m_lp.calcCoefficients( m_ac, Nc, m_extrap, condition ) != 0 ) + { + return internal::mxlib_error_report( error_t::liberr, "linearPredictor::calcCoefficients failed" ); + } + return error_t::noerror; } /// Worker function for regularizing the PSD for coefficient calculation. @@ -128,16 +159,19 @@ struct clAOLinearPredictor * * On subsequent calls, when min_var and min_sc are passed back in * loop over scale factors from min_sc-precision to max_sc in steps of + * + * \returns `error_t::noerror` on success, otherwise the coefficient-calculation error. */ template - int _regularizeCoefficients( realT &min_var, ///< [in.out] the minimum variance found. Set to 0 on initial call - realT &min_sc, ///< [in.out] the scale factor at the minimum variance. - realT precision, ///< [in] the step-size for the scale factor - realT max_sc, ///< [in] the maximum scale factor to test - clGainOpt &go_lp, ///< [in] the gain optimization object - std::vector &PSDt, ///< [in] the turbulence PSD - std::vector &PSDn, ///< [in] the WFS noise PSD - int Nc ///< [in] the number of coefficients + mx::error_t + _regularizeCoefficients( realT &min_var, /**< [in,out] the minimum variance found; set to 0 on initial call */ + realT &min_sc, /**< [in,out] the scale factor at the minimum variance */ + realT precision, /**< [in] the step size for the scale factor */ + realT max_sc, /**< [in] the maximum scale factor to test */ + clGainOpt &go_lp, /**< [in] the gain optimization object */ + std::vector &PSDt, /**< [in] the turbulence PSD */ + std::vector &PSDn, /**< [in] the WFS noise PSD */ + int Nc /**< [in] the number of coefficients */ ) { CLAOLP_BREADCRUMB; @@ -148,7 +182,7 @@ struct clAOLinearPredictor realT sc0; - if( min_var == 0 ) //first call + if( min_var == 0 ) // first call { sc0 = min_sc; min_var = std::numeric_limits::max(); @@ -165,13 +199,15 @@ struct clAOLinearPredictor CLAOLP_BREADCRUMB; // Test from sc0 to max_sc in steps of precision - //for( realT sc = sc0; sc <= max_sc; sc += precision ) + // for( realT sc = sc0; sc <= max_sc; sc += precision ) for( realT sc = max_sc; sc >= sc0; sc -= precision ) { CLAOLP_BREADCRUMB; - int rv = calcCoefficients( PSDt, PSDn, psdReg * pow( 10, -sc / 10 ), Nc ); - if( rv < 0 ) + ++m_regularizationReport.evaluations; + error_t rv = calcCoefficients( PSDt, PSDn, psdReg * pow( 10, -sc / 10 ), Nc ); + if( rv != error_t::noerror ) { + m_regularizationReport.status = regularizationStatus::calculationFailure; return rv; } @@ -194,7 +230,7 @@ struct clAOLinearPredictor if( telem ) { - m_regResults.push_back({sc, gopt_lp, gmax_lp, var_lp}); + m_regResults.push_back( { sc, gopt_lp, gmax_lp, var_lp } ); } CLAOLP_BREADCRUMB; @@ -207,14 +243,14 @@ struct clAOLinearPredictor // A jump by a factor of 10 indicates the wall if( var_lp > 10 * min_var ) { - return 0; + return error_t::noerror; } CLAOLP_BREADCRUMB; } CLAOLP_BREADCRUMB; - return 0; + return error_t::noerror; } /// Regularize the PSD and calculate the associated LP coefficients. @@ -222,29 +258,45 @@ struct clAOLinearPredictor * residual PSD. * * \tparam telem if true then the results are collected in m_regResults + * + * \returns `error_t::noerror` for a converged or boundary-limited search, `error_t::invalidconfig` for invalid + * controls, `error_t::timeout` on iteration exhaustion, or the coefficient-calculation error. */ template - int regularizeCoefficients( realT &gmax_lp, ///< [out] the maximum gain calculated for the regularized PSD - realT &gopt_lp, ///< [out] the optimum gain calculated for the regularized PSD - realT &var_lp, ///< [out] the variance at the optimum gain. - realT &min_sc, ///< [out] the optimum regularization scale factor - clGainOpt &go_lp, ///< [in] the gain optimization object - std::vector &PSDt, ///< [in] the turbulence PSD - std::vector &PSDn, ///< [in] the WFS noise PSD - int Nc ///< [in] the number of coefficients + mx::error_t + regularizeCoefficients( realT &gmax_lp, /**< [out] the maximum gain calculated for the regularized PSD */ + realT &gopt_lp, /**< [out] the optimum gain calculated for the regularized PSD */ + realT &var_lp, /**< [out] the variance at the optimum gain */ + realT &min_sc, /**< [out] the optimum regularization scale factor */ + clGainOpt &go_lp, /**< [in] the gain optimization object */ + std::vector &PSDt, /**< [in] the turbulence PSD */ + std::vector &PSDn, /**< [in] the WFS noise PSD */ + int Nc /**< [in] the number of coefficients */ ) { - CLAOLP_BREADCRUMB; + m_regularizationReport = {}; + + const realT intervalWidth = m_max_sc0 - m_min_sc0; + if( !math::isFinite( m_min_sc0 ) || !math::isFinite( m_max_sc0 ) || !math::isFinite( m_precision0 ) || + !math::isFinite( m_minPrecision ) || !math::isFinite( m_dPrecision ) || m_minPrecision <= 0 || + m_precision0 <= m_minPrecision || intervalWidth <= 0 || m_precision0 > intervalWidth || m_dPrecision <= 1 || + m_maxIts <= 0 ) + { + m_regularizationReport.status = regularizationStatus::invalidControls; + return internal::mxlib_error_report( error_t::invalidconfig, + "invalid linear-predictor regularization search controls" ); + } + realT min_var = m_min_var0; min_sc = m_min_sc0; realT precision = m_precision0; realT max_sc = m_max_sc0; - if(telem) + if( telem ) { - m_regResults.reserve(m_maxIts * 50); + m_regResults.reserve( m_maxIts * 50 ); } CLAOLP_BREADCRUMB; @@ -252,8 +304,11 @@ struct clAOLinearPredictor while( precision > m_minPrecision && its < m_maxIts ) { CLAOLP_BREADCRUMB; - int rv = _regularizeCoefficients( min_var, min_sc, precision, max_sc, go_lp, PSDt, PSDn, Nc ); - if( rv < 0) + const bool firstIteration = its == 0; + error_t rv = _regularizeCoefficients( min_var, min_sc, precision, max_sc, go_lp, PSDt, PSDn, Nc ); + ++its; + m_regularizationReport.iterations = its; + if( rv != error_t::noerror ) { return rv; } @@ -261,15 +316,14 @@ struct clAOLinearPredictor CLAOLP_BREADCRUMB; if( min_sc == max_sc ) { - if( its == 0 ) + if( firstIteration ) { min_sc -= precision; max_sc = 200; } else { - // std::cerr << "Error in regularizeCoefficients.\n"; - // return -1; + m_regularizationReport.status = regularizationStatus::boundaryLimited; break; } } @@ -278,15 +332,27 @@ struct clAOLinearPredictor max_sc = min_sc + precision; precision /= m_dPrecision; } + } - ++its; + if( precision > m_minPrecision && its >= m_maxIts && + m_regularizationReport.status != regularizationStatus::boundaryLimited ) + { + m_regularizationReport.status = regularizationStatus::iterationLimit; + return internal::mxlib_error_report( error_t::timeout, + "linear-predictor regularization reached its iteration limit" ); + } + + if( m_regularizationReport.status != regularizationStatus::boundaryLimited ) + { + m_regularizationReport.status = regularizationStatus::converged; } CLAOLP_BREADCRUMB; // Now record final values - int rv = calcCoefficients( PSDt, PSDn, PSDt[0] * pow( 10, -min_sc / 10 ), Nc ); - if( rv < 0 ) + error_t rv = calcCoefficients( PSDt, PSDn, PSDt[0] * pow( 10, -min_sc / 10 ), Nc ); + if( rv != error_t::noerror ) { + m_regularizationReport.status = regularizationStatus::calculationFailure; return rv; } @@ -300,25 +366,27 @@ struct clAOLinearPredictor gopt_lp = go_lp.optGainOpenLoop( var_lp, PSDt, PSDn, gmax_lp, false ); CLAOLP_BREADCRUMB; - return 0; + return error_t::noerror; } /// Regularize the PSD and calculate the associated LP coefficients. /** The PSD is regularized by adding a constant to it. This constant is found by minimizing the variance of the * residual PSD. * - * \tparam printout if true then the results are printed to stdout as they are calculated. + * \tparam printout if true then per-scale results are collected in m_regResults. + * + * \returns `error_t::noerror` on success, otherwise the regularization error. */ template - int optimizeNc( realT &gmax_lp, ///< [out] the maximum gain calculated for the regularized PSD - realT &gopt_lp, ///< [out] the optimum gain calculated for the regularized PSD - int &Nc, - realT &var_lp, ///< [out] the variance at the optimum gain. - clGainOpt &go_lp, ///< [in] the gain optimization object - std::vector &PSDt, ///< [in] the turbulence PSD - std::vector &PSDn, ///< [in] the WFS noise PSD - int minNc, ///< [in] the number of coefficients - int maxNc ) + mx::error_t optimizeNc( realT &gmax_lp, /**< [out] maximum gain for the selected predictor */ + realT &gopt_lp, /**< [out] optimum gain for the selected predictor */ + int &Nc, /**< [out] selected number of coefficients */ + realT &var_lp, /**< [out] variance at the optimum gain */ + clGainOpt &go_lp, /**< [in] the gain optimization object */ + std::vector &PSDt, /**< [in] the turbulence PSD */ + std::vector &PSDn, /**< [in] the WFS noise PSD */ + int minNc, /**< [in] minimum number of coefficients */ + int maxNc /**< [in] maximum number of coefficients */ ) { realT minVar = std::numeric_limits::max(); @@ -327,7 +395,12 @@ struct clAOLinearPredictor realT _gmax_lp; realT _gopt_lp; realT _var_lp; - regularizeCoefficients( _gmax_lp, _gopt_lp, _var_lp, go_lp, PSDt, PSDn, n ); + realT min_sc; + error_t rv = regularizeCoefficients( _gmax_lp, _gopt_lp, _var_lp, min_sc, go_lp, PSDt, PSDn, n ); + if( rv != error_t::noerror ) + { + return rv; + } if( _var_lp < minVar ) { @@ -340,7 +413,7 @@ struct clAOLinearPredictor } } - return 0; + return error_t::noerror; } }; diff --git a/include/ao/analysis/fourierTemporalPSD.hpp b/include/ao/analysis/fourierTemporalPSD.hpp index de0aea1fd..048bace8c 100644 --- a/include/ao/analysis/fourierTemporalPSD.hpp +++ b/include/ao/analysis/fourierTemporalPSD.hpp @@ -1139,20 +1139,20 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( doLP ) { realT min_sc; - int rv = tflp.regularizeCoefficients( gmax_lp, - gopt_lp, - var_lp, - min_sc, - go_lp, - tPSDpPOL, - tPSDn, - lpNc ); - - if( rv < 0 ) + error_t rv = tflp.regularizeCoefficients( gmax_lp, + gopt_lp, + var_lp, + min_sc, + go_lp, + tPSDpPOL, + tPSDn, + lpNc ); + + if( rv != error_t::noerror ) { std::cerr << "fourierTemporalPSD::analyzePSDGrid: regularizeCoefficients returned error "; - std::cerr << rv << ' '; + std::cerr << errorName( rv ) << ' '; std::cerr << __FILE__ << ' ' << __LINE__ << '\n'; } diff --git a/include/improc/imageUtils.hpp b/include/improc/imageUtils.hpp index bda8217ff..609940324 100644 --- a/include/improc/imageUtils.hpp +++ b/include/improc/imageUtils.hpp @@ -27,9 +27,10 @@ #ifndef improc_imageUtils_hpp #define improc_imageUtils_hpp -#include #include +#include "../math/floatUtils.hpp" + #include "imageTransforms.hpp" namespace mx @@ -77,13 +78,14 @@ constexpr T invalidNumber() return -3e38; } -/// Check if the number is nan, using several different methods -/** +/// Check whether a value represents an invalid image pixel. +/** Detects the mxlib invalid-number sentinel as well as NaN and positive or negative infinity. + * + * \returns true if value is invalid, otherwise false. */ -inline bool IsNan( float value ) +inline bool isInvalidPixel( float value /**< [in] value to test */ ) { - return ( ( ( ( *(uint32_t *)&value ) & 0x7fffffff ) > 0x7f800000 ) || ( value == invalidNumber() ) || - !std::isfinite( value ) ); + return value == invalidNumber() || !math::isFinite( value ); } /// Reflect pixel coordinates across the given center pixel. @@ -121,7 +123,7 @@ void zeroNaNs( imageT &im, ///< [in.out] image which will have any NaN pixels se { for( int r = 0; r < im.rows(); ++r ) { - if( IsNan( im( r, c ) ) ) + if( isInvalidPixel( im( r, c ) ) ) { im( r, c ) = val; } @@ -162,7 +164,7 @@ void zeroNaNCube( cubeT &imc, /**< [in.out] cube which will have any NaN pix { for( int r = 0; r < imc.rows(); ++r ) { - if( IsNan( imc.image( p )( r, c ) ) ) + if( isInvalidPixel( imc.image( p )( r, c ) ) ) { imc.image( p )( r, c ) = 0; if( mask ) @@ -333,11 +335,11 @@ imageMedian( const imageT &mat, /**< [in] the image */ template typename imageT::Scalar imageMedian( const imageT &mat, /**< [in] the image to take the median of*/ - std::vector *work = 0 /**< [in] [optional] working memory can - be retained and re-passed.*/ + std::vector *work = 0 /**< [in] [optional] working memory + can be retained and re-passed.*/ ) { - return imageMedian( mat, static_cast *>(nullptr), work ); + return imageMedian( mat, static_cast *>( nullptr ), work ); } /// Calculate the center of light of an image @@ -528,7 +530,6 @@ void removeCols( eigenT &out, const eigenTin &in, int st, int w ) out.topRightCorner( in.rows(), in.cols() - ( st + w ) ) = in.topRightCorner( in.rows(), in.cols() - ( st + w ) ); } - /** \ingroup image_utils *@{ */ @@ -581,8 +582,6 @@ void *imcpy_flipUDLR( void *dest, ///< [out] the address of the first pixel i size_t szof ///< [in] the size in bytes of a one pixel ); - - } // namespace improc } // namespace mx diff --git a/include/math/CMakeLists.txt b/include/math/CMakeLists.txt index 91606f9c9..878f2ad68 100644 --- a/include/math/CMakeLists.txt +++ b/include/math/CMakeLists.txt @@ -6,6 +6,7 @@ add_subdirectory(plot) set(OBJLIB_INCLUDES ${OBJLIB_INCLUDES} include/math/constants.hpp include/math/eigenLapack.hpp + include/math/floatUtils.hpp include/math/geo.hpp include/math/gslInterpolation.hpp include/math/gslInterpolator.hpp diff --git a/include/math/floatUtils.hpp b/include/math/floatUtils.hpp new file mode 100644 index 000000000..feca7051d --- /dev/null +++ b/include/math/floatUtils.hpp @@ -0,0 +1,95 @@ +/** \file floatUtils.hpp + * \author Jared R. Males + * \brief Floating-point classification utilities that remain reliable under fast-math optimization. + * \ingroup gen_math_files + */ + +#ifndef math_floatUtils_hpp +#define math_floatUtils_hpp + +#include +#include +#include +#include + +namespace mx +{ +namespace math +{ + +namespace floatUtils_detail +{ + +/// Convert an extended floating-point value to a classifiable double without overflowing finite values. +template +double normalizedDouble( realT value /**< [in] floating-point value to normalize */ ) +{ + return static_cast( value / std::numeric_limits::max() ); +} + +} // namespace floatUtils_detail + +/// Test whether a floating-point value is NaN, including under finite-math-only optimization. +/** + * \returns true if value is a quiet or signaling NaN, otherwise false. + * + * \ingroup gen_math + */ +template +bool isNan( realT value /**< [in] floating-point value to test */ ) +{ + static_assert( std::is_floating_point_v, "isNan requires a floating-point type" ); + + if constexpr( std::numeric_limits::is_iec559 && sizeof( realT ) == sizeof( std::uint32_t ) ) + { + constexpr std::uint32_t exponentMask = 0x7f800000U; + constexpr std::uint32_t mantissaMask = 0x007fffffU; + const std::uint32_t bits = std::bit_cast( value ); + return ( bits & exponentMask ) == exponentMask && ( bits & mantissaMask ) != 0; + } + else if constexpr( std::numeric_limits::is_iec559 && sizeof( realT ) == sizeof( std::uint64_t ) ) + { + constexpr std::uint64_t exponentMask = 0x7ff0000000000000ULL; + constexpr std::uint64_t mantissaMask = 0x000fffffffffffffULL; + const std::uint64_t bits = std::bit_cast( value ); + return ( bits & exponentMask ) == exponentMask && ( bits & mantissaMask ) != 0; + } + else + { + const double normalized = floatUtils_detail::normalizedDouble( value ); + return isNan( normalized ); + } +} + +/// Test whether a floating-point value is finite, including under finite-math-only optimization. +/** + * \returns true if value is neither infinite nor NaN, otherwise false. + * + * \ingroup gen_math + */ +template +bool isFinite( realT value /**< [in] floating-point value to test */ ) +{ + static_assert( std::is_floating_point_v, "isFinite requires a floating-point type" ); + + if constexpr( std::numeric_limits::is_iec559 && sizeof( realT ) == sizeof( std::uint32_t ) ) + { + constexpr std::uint32_t exponentMask = 0x7f800000U; + return ( std::bit_cast( value ) & exponentMask ) != exponentMask; + } + else if constexpr( std::numeric_limits::is_iec559 && sizeof( realT ) == sizeof( std::uint64_t ) ) + { + constexpr std::uint64_t exponentMask = 0x7ff0000000000000ULL; + return ( std::bit_cast( value ) & exponentMask ) != exponentMask; + } + else + { + const double normalized = floatUtils_detail::normalizedDouble( value ); + return isFinite( normalized ); + } +} + +} // namespace math +} // namespace mx + +#endif // math_floatUtils_hpp diff --git a/include/math/math.hpp b/include/math/math.hpp index 7e9f49d74..ccf083cec 100644 --- a/include/math/math.hpp +++ b/include/math/math.hpp @@ -53,6 +53,7 @@ #include "plot/gnuPlot.hpp" #include "constants.hpp" #include "eigenLapack.hpp" +#include "floatUtils.hpp" #include "geo.hpp" #include "gslInterpolation.hpp" #include "gslInterpolator.hpp" diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fc6ba464b..ede0af842 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,7 @@ set(MXLIB_TEST_MAIN_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/testMain.cpp) set(MXLIB_TEST_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/include/ao/analysis/aoAtmosphere_test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ao/analysis/clAOLinearPredictor_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ao/analysis/aoSystem_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ao/analysis/aoPSDs_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/astro/astroDynamics_test.cpp @@ -13,6 +14,7 @@ set(MXLIB_TEST_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/include/ioutils/fits/fitsHeaderCard_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ioutils/fits/fitsFile_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/math/geo_test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/math/floatUtils_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/math/func/moffat_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/math/templateBLAS_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/include/math/templateLapack_test.cpp diff --git a/tests/Makefile b/tests/Makefile index e068e49f8..83537fe55 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -6,6 +6,7 @@ INCLUDES += -I../include OBJS = testMain.o \ include/ao/analysis/aoAtmosphere_test.o \ + include/ao/analysis/clAOLinearPredictor_test.o \ include/ao/analysis/aoSystem_test.o \ include/ao/analysis/aoPSDs_test.o \ include/astro/astroDynamics_test.o \ @@ -13,6 +14,7 @@ OBJS = testMain.o \ include/ioutils/fits/fitsHeaderCard_test.o \ include/ioutils/fits/fitsFile_test.o \ include/math/geo_test.o \ + include/math/floatUtils_test.o \ include/math/func/moffat_test.o \ include/math/templateBLAS_test.o \ include/math/templateLapack_test.o \ diff --git a/tests/include/ao/analysis/clAOLinearPredictor_test.cpp b/tests/include/ao/analysis/clAOLinearPredictor_test.cpp new file mode 100644 index 000000000..08620b55b --- /dev/null +++ b/tests/include/ao/analysis/clAOLinearPredictor_test.cpp @@ -0,0 +1,209 @@ +/** \file clAOLinearPredictor_test.cpp + * \brief Tests of closed-loop linear-predictor regularization. + */ + +#include "../../../catch2/catch.hpp" + +#define MX_NO_ERROR_REPORTS + +#include "../../../../include/ao/analysis/clAOLinearPredictor.hpp" + +#include +#include +#include + +namespace +{ + +using predictorT = mx::AO::analysis::clAOLinearPredictor; +using optimizerT = mx::AO::analysis::clGainOpt; + +/// Construct a compact, valid PSD fixture for regularization tests. +void makeFixture( optimizerT &optimizer, /**< [out] configured gain optimizer */ + std::vector &disturbance, /**< [out] positive disturbance PSD */ + std::vector &noise /**< [out] nonnegative noise PSD */ ) +{ + constexpr std::size_t sampleCount = 64; + std::vector frequency( sampleCount ); + disturbance.resize( sampleCount ); + noise.assign( sampleCount, 1e-3 ); + + for( std::size_t index = 0; index < sampleCount; ++index ) + { + frequency[index] = 0.5 * static_cast( index + 1 ) / static_cast( sampleCount ); + disturbance[index] = 1.0 / ( 1.0 + 100.0 * frequency[index] * frequency[index] ); + } + + optimizer.f( frequency ); +} + +/// Run one regularization search with telemetry enabled. +mx::error_t runSearch( predictorT &predictor, /**< [in,out] predictor under test */ + optimizerT &optimizer, /**< [in,out] configured gain optimizer */ + std::vector &disturbance, /**< [in] disturbance PSD */ + std::vector &noise /**< [in] noise PSD */ ) +{ + double maximumGain = 0; + double optimalGain = 0; + double variance = 0; + double scale = 0; + return predictor + .regularizeCoefficients( maximumGain, optimalGain, variance, scale, optimizer, disturbance, noise, 4 ); +} + +} // namespace + +/// Verify that invalid linear-predictor regularization controls are rejected before evaluation. +/** Exercises validation of the regularization interval, spacing, refinement divisor, and iteration limit. */ +SCENARIO( "Linear-predictor regularization rejects invalid search controls", "[ao::analysis::clAOLinearPredictor]" ) +{ + // clang-format off +#ifdef __DOXY_ONLY__ + mx::AO::analysis::clAOLinearPredictor::regularizeCoefficients(); +#endif + // clang-format on + + optimizerT optimizer( 1.0, 1.5 ); + std::vector disturbance; + std::vector noise; + makeFixture( optimizer, disturbance, noise ); + + GIVEN( "a fresh predictor" ) + { + predictorT predictor; + + WHEN( "the initial precision is below the minimum" ) + { + predictor.m_precision0 = 0.5 * predictor.m_minPrecision; + REQUIRE( runSearch( predictor, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + } + + WHEN( "the initial precision equals the minimum" ) + { + predictor.m_precision0 = predictor.m_minPrecision; + REQUIRE( runSearch( predictor, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + } + + WHEN( "the initial precision is zero, negative, NaN, or infinite" ) + { + const std::vector invalidValues{ 0, + -1, + std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity() }; + for( const double value : invalidValues ) + { + predictorT invalidPredictor; + invalidPredictor.m_precision0 = value; + REQUIRE( runSearch( invalidPredictor, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + REQUIRE( invalidPredictor.m_regularizationReport.status == + predictorT::regularizationStatus::invalidControls ); + REQUIRE( invalidPredictor.m_regularizationReport.evaluations == 0 ); + REQUIRE( invalidPredictor.m_regResults.empty() ); + } + } + + WHEN( "the initial precision is wider than the initial interval" ) + { + predictor.m_precision0 = predictor.m_max_sc0 - predictor.m_min_sc0 + 1; + REQUIRE( runSearch( predictor, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + } + + WHEN( "the scale interval, refinement divisor, or iteration limit is invalid" ) + { + predictorT reversedInterval; + reversedInterval.m_max_sc0 = reversedInterval.m_min_sc0; + REQUIRE( runSearch( reversedInterval, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + + predictorT invalidDivisor; + invalidDivisor.m_dPrecision = 1; + REQUIRE( runSearch( invalidDivisor, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + + predictorT invalidIterations; + invalidIterations.m_maxIts = 0; + REQUIRE( runSearch( invalidIterations, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + } + + WHEN( "another floating-point search control is nonfinite or nonpositive" ) + { + predictorT nonfiniteMinimumScale; + nonfiniteMinimumScale.m_min_sc0 = std::numeric_limits::quiet_NaN(); + REQUIRE( runSearch( nonfiniteMinimumScale, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + + predictorT nonfiniteMaximumScale; + nonfiniteMaximumScale.m_max_sc0 = std::numeric_limits::infinity(); + REQUIRE( runSearch( nonfiniteMaximumScale, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + + predictorT zeroMinimumPrecision; + zeroMinimumPrecision.m_minPrecision = 0; + REQUIRE( runSearch( zeroMinimumPrecision, optimizer, disturbance, noise ) == mx::error_t::invalidconfig ); + + predictorT nonfiniteRefinementDivisor; + nonfiniteRefinementDivisor.m_dPrecision = std::numeric_limits::infinity(); + REQUIRE( runSearch( nonfiniteRefinementDivisor, optimizer, disturbance, noise ) == + mx::error_t::invalidconfig ); + } + } +} + +/// Verify that linear-predictor regularization reports how its search terminated. +/** Exercises completed, boundary-limited, and iteration-limited regularization searches. */ +SCENARIO( "Linear-predictor regularization reports search termination", "[ao::analysis::clAOLinearPredictor]" ) +{ + // clang-format off +#ifdef __DOXY_ONLY__ + mx::AO::analysis::clAOLinearPredictor::regularizeCoefficients(); +#endif + // clang-format on + + optimizerT optimizer( 1.0, 1.5 ); + std::vector disturbance; + std::vector noise; + makeFixture( optimizer, disturbance, noise ); + + GIVEN( "a precision immediately above the minimum and one allowed iteration" ) + { + predictorT predictor; + predictor.m_min_sc0 = 10; + predictor.m_max_sc0 = 10.0025; + predictor.m_precision0 = std::nextafter( predictor.m_minPrecision, std::numeric_limits::infinity() ); + predictor.m_maxIts = 1; + + const mx::error_t result = runSearch( predictor, optimizer, disturbance, noise ); + + REQUIRE( result != mx::error_t::invalidconfig ); + REQUIRE( predictor.m_regularizationReport.status != predictorT::regularizationStatus::invalidControls ); + REQUIRE( predictor.m_regularizationReport.iterations == 1 ); + REQUIRE( predictor.m_regularizationReport.evaluations > 0 ); + REQUIRE_FALSE( predictor.m_regResults.empty() ); + } + + GIVEN( "a valid search forced to exhaust its iteration limit" ) + { + predictorT predictor; + predictor.m_min_sc0 = 10; + predictor.m_max_sc0 = 12; + predictor.m_precision0 = 1; + predictor.m_minPrecision = 1e-9; + predictor.m_maxIts = 1; + + REQUIRE( runSearch( predictor, optimizer, disturbance, noise ) == mx::error_t::timeout ); + REQUIRE( predictor.m_regularizationReport.status == predictorT::regularizationStatus::iterationLimit ); + REQUIRE( predictor.m_regularizationReport.iterations == 1 ); + REQUIRE( predictor.m_regularizationReport.evaluations > 0 ); + } + + GIVEN( "a valid search whose initial precision equals the interval width" ) + { + predictorT predictor; + predictor.m_min_sc0 = 10; + predictor.m_max_sc0 = 12; + predictor.m_precision0 = 2; + predictor.m_minPrecision = 0.1; + + REQUIRE( runSearch( predictor, optimizer, disturbance, noise ) == mx::error_t::noerror ); + REQUIRE( ( predictor.m_regularizationReport.status == predictorT::regularizationStatus::converged || + predictor.m_regularizationReport.status == predictorT::regularizationStatus::boundaryLimited ) ); + REQUIRE( predictor.m_regularizationReport.iterations > 0 ); + REQUIRE( predictor.m_regularizationReport.evaluations > 0 ); + } +} diff --git a/tests/include/improc/imageUtils_test.cpp b/tests/include/improc/imageUtils_test.cpp index 134f35f02..88b5ea2aa 100644 --- a/tests/include/improc/imageUtils_test.cpp +++ b/tests/include/improc/imageUtils_test.cpp @@ -1,16 +1,30 @@ /** \file imageUtils_test.cpp + * \brief Tests of image-processing utilities. */ #include "../../catch2/catch.hpp" -#include #include +#include +#include + #define MX_NO_ERROR_REPORTS #include "../../../include/math/func/gaussian.hpp" #include "../../../include/improc/imageUtils.hpp" #include "../../../include/improc/eigenCube.hpp" +/// Verify image invalid-pixel classification for the sentinel and nonfinite values. +/** Preserves the established invalid-pixel contract. */ +TEST_CASE( "Image invalid-pixel detection handles sentinel and nonfinite values", "[improc::isInvalidPixel]" ) +{ + REQUIRE_FALSE( mx::improc::isInvalidPixel( 0.0F ) ); + REQUIRE( mx::improc::isInvalidPixel( std::numeric_limits::quiet_NaN() ) ); + REQUIRE( mx::improc::isInvalidPixel( std::numeric_limits::infinity() ) ); + REQUIRE( mx::improc::isInvalidPixel( -std::numeric_limits::infinity() ) ); + REQUIRE( mx::improc::isInvalidPixel( mx::improc::invalidNumber() ) ); +} + /** Scenario: centroiding Gaussians with center of light * * Verify center of light calculation diff --git a/tests/include/math/floatUtils_test.cpp b/tests/include/math/floatUtils_test.cpp new file mode 100644 index 000000000..fa4d53ae0 --- /dev/null +++ b/tests/include/math/floatUtils_test.cpp @@ -0,0 +1,31 @@ +/** \file floatUtils_test.cpp + * \brief Tests of floating-point classification utilities. + */ + +#include "../../catch2/catch.hpp" + +#include "../../../include/math/floatUtils.hpp" + +#include + +/// Verify floating-point classification for finite and nonfinite values. +/** Exercises the NaN and finite-value classifiers for float, double, and long double. + * The same test applies to standard and fast-math builds. + */ +TEST_CASE( "Floating-point classification handles finite and nonfinite values", "[math::floatUtils]" ) +{ + REQUIRE( mx::math::isFinite( std::numeric_limits::max() ) ); + REQUIRE_FALSE( mx::math::isFinite( std::numeric_limits::infinity() ) ); + REQUIRE( mx::math::isNan( std::numeric_limits::quiet_NaN() ) ); + REQUIRE_FALSE( mx::math::isNan( std::numeric_limits::infinity() ) ); + + REQUIRE( mx::math::isFinite( std::numeric_limits::max() ) ); + REQUIRE_FALSE( mx::math::isFinite( std::numeric_limits::quiet_NaN() ) ); + REQUIRE( mx::math::isNan( std::numeric_limits::quiet_NaN() ) ); + REQUIRE_FALSE( mx::math::isNan( std::numeric_limits::infinity() ) ); + + REQUIRE( mx::math::isFinite( std::numeric_limits::max() ) ); + REQUIRE_FALSE( mx::math::isFinite( std::numeric_limits::infinity() ) ); + REQUIRE( mx::math::isNan( std::numeric_limits::quiet_NaN() ) ); + REQUIRE_FALSE( mx::math::isNan( std::numeric_limits::infinity() ) ); +} From 1e8c4bb203db0a27fcd9f6ec39fc1d6fdaa3d5de Mon Sep 17 00:00:00 2001 From: Jared Males Date: Sun, 2 Aug 2026 10:54:07 -0700 Subject: [PATCH 2/2] Document automatic Doxygen test references --- AGENTS.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/AGENTS.md b/AGENTS.md index e249807fc..52a5dd4d1 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -59,6 +59,19 @@ Follow these code style and documentation rules exactly. - `This work was performed by GPT-5.3-Codex in response to the prompt: "...".` - Include the primary user prompt verbatim (or a faithful condensed version if it is extremely long). +12) Unit Test Documentation +- Add a brief Doxygen block immediately before every Catch2 `TEST_CASE` or `SCENARIO`. +- State the behavior being verified and identify the real production API under test. +- Let Doxygen discover real calls in the test body so the test appears in each production API's `Referenced by` list. +- Do not use `\test` or prose-only `\ref` commands to manufacture test-to-API links. + +13) Preserve Doxygen Links Through Test Harnesses +- Preserve Doxygen links to the real production APIs when test fixtures, wrappers, namespaces, macros, or private-access techniques prevent automatic symbol linking. +- Add explicit Doxygen-only code references to the production symbols inside the relevant test body when direct calls are otherwise hidden. +- Guard reference-only code with `#ifdef __DOXY_ONLY__` so it need not compile, and use raw calls or member references that Doxygen can add to the production symbol's `Referenced by` list. +- Hide harness-only helpers from generated documentation with `\cond` and `\endcond` when they would dominate or obscure production API links. +- Disable `clang-format` around non-compiling Doxygen-only reference blocks when necessary. + When you finish: - Summarize what changed. - List affected files.