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. diff --git a/include/ao/analysis/aoAtmosphere.hpp b/include/ao/analysis/aoAtmosphere.hpp index 0d45f890b..a551eaa0f 100644 --- a/include/ao/analysis/aoAtmosphere.hpp +++ b/include/ao/analysis/aoAtmosphere.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -20,6 +21,7 @@ #include "aoConstants.hpp" #include "../../math/constants.hpp" +#include "../../math/floatUtils.hpp" #include "../../app/appConfigurator.hpp" @@ -49,51 +51,52 @@ class aoAtmosphere /// Constructor aoAtmosphere(); + /// Validate the complete atmosphere configuration for calculation. + /** Incremental setters may temporarily leave layer vectors incomplete. Call this at configuration and calculation + * boundaries before using indexed or derived atmosphere values. + * + * \returns `error_t::noerror` when all scalar and layer invariants are satisfied, or a typed configuration error. + */ + error_t validate() const; + protected: - realT m_r_0; ///< Fried's parameter, in m + realT m_r_0{ 0 }; ///< Fried's parameter, in m - realT m_lam_0{ 0.5e-6 }; ///< Wavelength of Fried's parameter, in m + realT m_lam_0{ 0.5e-6 }; ///< Wavelength of Fried's parameter, in m - std::vector m_layer_Cn2; ///< Vector of layer strengths. + std::vector m_layer_Cn2; ///< Vector of layer strengths. - std::vector m_L_0; ///< The outer scale, in m + std::vector m_L_0; ///< The outer scale, in m - std::vector m_l_0; ///< The inner scale of each layer, in m + std::vector m_l_0; ///< The inner scale of each layer, in m - bool m_nonKolmogorov{ false }; ///< Flag indicating if non-Kolmogorov PSD parameters are used. + bool m_nonKolmogorov{ false }; ///< Flag indicating if non-Kolmogorov PSD parameters are used. - std::vector m_beta{ 1 }; ///< The PSD normalization when in non-Kolmogorov mode. + std::vector m_beta{ 1 }; ///< The PSD normalization when in non-Kolmogorov mode. - std::vector m_alpha{ 0 }; ///< The PSD exponent when in non-Kolmogorov mode. + std::vector m_alpha{ 0 }; ///< The PSD exponent when in non-Kolmogorov mode. - std::vector m_beta_0{ 0 }; ///< The PSD constant when in non-Kolmogorov mode. + std::vector m_beta_0{ 0 }; ///< The PSD constant when in non-Kolmogorov mode. - std::vector m_layer_z; ///< Vector of layer heights, in m, above the observatory. + std::vector m_layer_z; ///< Vector of layer heights, in m, above the observatory. - realT m_h_obs{ 0 }; ///< Height of the observatory above sea level, in m. + realT m_h_obs{ 0 }; ///< Height of the observatory above sea level, in m. - realT m_H{ 8000 }; ///< The atmospheric scale height, in m. + realT m_H{ 8000 }; ///< The atmospheric scale height, in m. std::vector m_layer_v_wind; ///< Vector of layer wind speeds, in m/s. - std::vector m_layer_dir; ///< Vector of layer wind directions, in radians. + std::vector m_layer_dir; ///< Vector of layer wind directions, in radians. - bool m_v_wind_updated{ false }; ///< whether or not m_v_wind has been updated after changes + bool m_v_wind_updated{ false }; ///< whether or not m_v_wind has been updated after changes - realT m_v_wind; ///< \f$ C_n^2 \f$ averaged windspeed + realT m_v_wind{ 0 }; ///< \f$ C_n^2 \f$ averaged windspeed - realT m_dir_wind; ///< \f$ C_n^2 \f$ averaged direction + realT m_dir_wind{ 0 }; ///< \f$ C_n^2 \f$ averaged direction - bool m_z_mean_updated{ false }; ///< whether or not m_z_mean has been updated after changes + bool m_z_mean_updated{ false }; ///< whether or not m_z_mean has been updated after changes - realT m_z_mean; ///< \f$ C_n^2 \f$ averaged layer height - - /// Checks if layer vectors have consistent length. - /** - * \returns 0 if all layer vectors are the same length - * \returns -1 if not, and prints an error. - */ - int checkLayers(); + realT m_z_mean{ 0 }; ///< \f$ C_n^2 \f$ averaged layer height public: /** \name PSD Parameters @@ -150,8 +153,8 @@ class aoAtmosphere * \f$ \sum_n C_n^2 = 1 \f$. * */ - void layer_Cn2( const std::vector &cn2, ///< [in] is a vector containing the layer strengths - const realT l0 = 0 ///< [in] [optional] if l0 > 0, then r_0 is set from the layer strengths. + error_t layer_Cn2( const std::vector &cn2, ///< [in] vector containing the layer strengths + const realT l0 = 0 ///< [in] reference wavelength; if positive, also calculate r_0 ); /// Get the value of the outer scale for a single layer. @@ -598,19 +601,15 @@ class aoAtmosphere /** @{ */ - /// Setup the configurator to configure this class - /** - * Tests: - * - Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]" - */ + /// Setup the configurator to configure this class. void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); - /// Load the configuration of this class from a configurator - /** - * Tests: - * - Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]" + /// Load the configuration of this class from a configurator. + /** The complete atmosphere is validated before derived rescalings are applied and again before returning. + * + * \returns `error_t::noerror` for a valid loaded atmosphere, or a typed configuration error. */ - void loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); + error_t loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); /// @} }; @@ -621,41 +620,63 @@ aoAtmosphere::aoAtmosphere() } template -int aoAtmosphere::checkLayers() +error_t aoAtmosphere::validate() const { - size_t n = m_L_0.size(); + const size_t layerCount = m_layer_Cn2.size(); + if( layerCount == 0 ) + { + return internal::mxlib_error_report( error_t::sizeerr, "atmosphere must contain at least one layer" ); + } - if( m_l_0.size() != n ) + if( m_L_0.size() != layerCount || m_l_0.size() != layerCount || m_layer_z.size() != layerCount || + m_layer_v_wind.size() != layerCount || m_layer_dir.size() != layerCount ) { - internal::mxlib_error_report(error_t::sizeerr,"mismatched layer numbers (inner scale vs. outer scale)"); - return -1; + return internal::mxlib_error_report( error_t::sizeerr, "atmosphere layer-vector sizes do not match" ); } - if( m_layer_z.size() != n ) + if( m_nonKolmogorov && + ( m_alpha.size() != layerCount || m_beta.size() != layerCount || m_beta_0.size() != layerCount ) ) { - internal::mxlib_error_report(error_t::sizeerr,"mismatched layer numbers (layer_z vs. outer scale)"); - return -1; + return internal::mxlib_error_report( error_t::sizeerr, "non-Kolmogorov atmosphere vector sizes do not match" ); } - if( m_layer_Cn2.size() != n ) + if( !math::isFinite( m_h_obs ) || m_h_obs < 0 || !math::isFinite( m_H ) || m_H <= 0 || + ( !m_nonKolmogorov && + ( !math::isFinite( m_r_0 ) || m_r_0 <= 0 || !math::isFinite( m_lam_0 ) || m_lam_0 <= 0 ) ) ) { - internal::mxlib_error_report(error_t::sizeerr,"mismatched layer numbers (layer_Cn2 vs. outer scale)" ); - return -1; + return internal::mxlib_error_report( error_t::invalidconfig, "atmosphere scalar parameters are invalid" ); } - if( m_layer_dir.size() != n ) + realT totalStrength = 0; + bool hasPositiveStrength = false; + for( size_t index = 0; index < layerCount; ++index ) { - internal::mxlib_error_report(error_t::sizeerr,"mismatched layer numbers (layer_dir vs. outer scale)"); - return -1; + if( !math::isFinite( m_layer_Cn2[index] ) || m_layer_Cn2[index] < 0 || !math::isFinite( m_L_0[index] ) || + !math::isFinite( m_l_0[index] ) || m_l_0[index] < 0 || !math::isFinite( m_layer_z[index] ) || + m_layer_z[index] < 0 || !math::isFinite( m_layer_v_wind[index] ) || m_layer_v_wind[index] < 0 || + !math::isFinite( m_layer_dir[index] ) ) + { + return internal::mxlib_error_report( error_t::invalidconfig, "atmosphere layer parameters are invalid" ); + } + + totalStrength += m_layer_Cn2[index]; + hasPositiveStrength = hasPositiveStrength || m_layer_Cn2[index] > 0; + + if( m_nonKolmogorov && ( !math::isFinite( m_alpha[index] ) || !math::isFinite( m_beta[index] ) || + m_beta[index] <= 0 || !math::isFinite( m_beta_0[index] ) || m_beta_0[index] < 0 ) ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "non-Kolmogorov atmosphere parameters are invalid" ); + } } - if( m_layer_v_wind.size() != n ) + if( !hasPositiveStrength || !math::isFinite( totalStrength ) || totalStrength <= 0 ) { - internal::mxlib_error_report(error_t::sizeerr,"mismatched layer numbers (layer_v_wind vs. outer scale)"); - return -1; + return internal::mxlib_error_report( error_t::invalidconfig, + "atmosphere must contain a finite positive layer-strength sum" ); } - return 0; + return error_t::noerror; } template @@ -704,28 +725,54 @@ std::vector aoAtmosphere::layer_Cn2() } template -void aoAtmosphere::layer_Cn2( const std::vector &cn2, const realT l0 ) +error_t aoAtmosphere::layer_Cn2( const std::vector &cn2, const realT l0 ) { - m_layer_Cn2 = cn2; + if( cn2.empty() || !math::isFinite( l0 ) || l0 < 0 ) + { + return internal::mxlib_error_report( error_t::invalidarg, + "layer strengths must be nonempty and reference wavelength nonnegative" ); + } realT layer_norm = 0; - - for( size_t i = 0; i < m_layer_Cn2.size(); ++i ) + for( size_t i = 0; i < cn2.size(); ++i ) { + if( !math::isFinite( cn2[i] ) || cn2[i] < 0 ) + { + return internal::mxlib_error_report( error_t::invalidarg, + "layer strengths must be finite and nonnegative" ); + } layer_norm += cn2[i]; } - for( size_t i = 0; i < m_layer_Cn2.size(); ++i ) - m_layer_Cn2[i] = m_layer_Cn2[i] / layer_norm; + if( !math::isFinite( layer_norm ) || layer_norm <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidarg, "layer strengths must have a finite positive sum" ); + } + + std::vector normalized = cn2; + for( size_t i = 0; i < normalized.size(); ++i ) + { + normalized[i] /= layer_norm; + } if( l0 > 0 ) { - m_r_0 = 1.0 / pow( layer_norm * 5.520e13, math::three_fifths() ); + const realT r0 = 1.0 / pow( layer_norm * 5.520e13, math::three_fifths() ); + if( !math::isFinite( r0 ) || r0 <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidarg, + "layer strengths produce an invalid Fried parameter" ); + } + m_r_0 = r0; m_lam_0 = l0; } + m_layer_Cn2 = std::move( normalized ); + m_v_wind_updated = false; m_z_mean_updated = false; + + return error_t::noerror; } template @@ -1211,6 +1258,8 @@ void aoAtmosphere::loadGuyon2005() layer_z( { 500, 1000, 2000, 4000, 8000, 16000 } ); layer_v_wind( { 10., 10., 10., 10., 10., 10. } ); layer_dir( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ); + L_0( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ); + l_0( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ); r_0( 0.2, 0.5e-6 ); @@ -1329,10 +1378,24 @@ void aoAtmosphere::setupConfig( app::appConfigurator &config ) false, "real", "The reference wavlength for r_0 [m]" ); - config.add( - "atm.L_0", "", "atm.L_0", argType::Required, "atm", "L_0", false, "vector", "Layer outer scales [m]" ); - config.add( - "atm.l_0", "", "atm.l_0", argType::Required, "atm", "l_0", false, "vector", "Layer inner scales [m]" ); + config.add( "atm.L_0", + "", + "atm.L_0", + argType::Required, + "atm", + "L_0", + false, + "vector", + "Layer outer scales [m]" ); + config.add( "atm.l_0", + "", + "atm.l_0", + argType::Required, + "atm", + "l_0", + false, + "vector", + "Layer inner scales [m]" ); config.add( "atm.layer_z", "", "atm.layer_z", @@ -1342,8 +1405,15 @@ void aoAtmosphere::setupConfig( app::appConfigurator &config ) false, "vector", "layer heights [m]" ); - config.add( - "atm.h_obs", "", "atm.h_obs", argType::Required, "atm", "h_obs", false, "real", "height of observatory [m]" ); + config.add( "atm.h_obs", + "", + "atm.h_obs", + argType::Required, + "atm", + "h_obs", + false, + "real", + "height of observatory [m]" ); config.add( "atm.H", "", "atm.H", argType::Required, "atm", "H", false, "real", "atmospheric scale heights [m]" ); config.add( "atm.layer_Cn2", "", @@ -1438,7 +1508,7 @@ void aoAtmosphere::setupConfig( app::appConfigurator &config ) } template -void aoAtmosphere::loadConfig( app::appConfigurator &config ) +error_t aoAtmosphere::loadConfig( app::appConfigurator &config ) { // Here "has side effecs" means that the set function does more than simply copy the value. @@ -1452,7 +1522,13 @@ void aoAtmosphere::loadConfig( app::appConfigurator &config ) std::vector lcn2 = m_layer_Cn2; config( lcn2, "atm.layer_Cn2" ); if( config.isSet( "atm.layer_Cn2" ) ) - layer_Cn2( lcn2 ); + { + const error_t strengthStatus = layer_Cn2( lcn2 ); + if( strengthStatus != error_t::noerror ) + { + return strengthStatus; + } + } realT r0 = r_0(); config( r0, "atm.r_0" ); @@ -1467,7 +1543,7 @@ void aoAtmosphere::loadConfig( app::appConfigurator &config ) std::vector layz = m_layer_z; config( layz, "atm.layer_z" ); // Do this no matter what to record source if( config.isSet( "atm.layer_z" ) ) - layer_z( layz ); // but only call this if changed + layer_z( layz ); // but only call this if changed config( m_h_obs, "atm.h_obs" ); config( m_H, "atm.H" ); @@ -1476,30 +1552,21 @@ void aoAtmosphere::loadConfig( app::appConfigurator &config ) std::vector lvw = m_layer_v_wind; config( lvw, "atm.layer_v_wind" ); // Do this no matter what to record source if( config.isSet( "atm.layer_v_wind" ) ) - layer_v_wind( lvw ); // but only call this if changed + layer_v_wind( lvw ); // but only call this if changed // Has side effects: std::vector ld = m_layer_dir; config( ld, "atm.layer_dir" ); // Do this no matter what to record source if( config.isSet( "atm.layer_dir" ) ) - layer_dir( ld ); // but only call this if changed + layer_dir( ld ); // but only call this if changed - realT vw = m_v_wind; + realT vw = 0; config( vw, "atm.v_wind" ); // Do this no matter what to record source - if( config.isSet( "atm.v_wind" ) ) - v_wind( vw ); // but only call this if changed - realT t0 = tau_0(); - config( t0, "atm.tau_0" ); // Do this no matter what to record source - if( config.isSet( "atm.tau_0" ) ) - { - std::cerr << "setting tau_0 " << t0 << "\n"; - tau_0( t0, m_lam_0 ); // but only call this if changed - } - realT zm = m_z_mean; + realT t0 = 0; + config( t0, "atm.tau_0" ); // Do this no matter what to record source + realT zm = 0; config( zm, "atm.z_mean" ); // Do this no matter what to record source - if( config.isSet( "atm.z_mean" ) ) - z_mean( zm ); // but only call this if changed config( m_nonKolmogorov, "atm.nonKolmogorov" ); @@ -1517,6 +1584,66 @@ void aoAtmosphere::loadConfig( app::appConfigurator &config ) config( b0, "atm.beta_0" ); if( config.isSet( "atm.beta_0" ) ) beta_0( b0 ); // this sets m_nonKolmogorov + + error_t status = validate(); + if( status != error_t::noerror ) + { + return status; + } + + if( config.isSet( "atm.v_wind" ) ) + { + if( !math::isFinite( vw ) || vw <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "configured mean wind speed must be finite and positive" ); + } + + if( v_wind() <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "a static atmosphere cannot be rescaled to positive mean wind" ); + } + v_wind( vw ); + } + + if( config.isSet( "atm.tau_0" ) ) + { + if( !math::isFinite( t0 ) || t0 <= 0 || !math::isFinite( m_r_0 ) || m_r_0 <= 0 || !math::isFinite( m_lam_0 ) || + m_lam_0 <= 0 || v_wind() <= 0 ) + { + return internal::mxlib_error_report( + error_t::invalidconfig, + "configured atmosphere time constant requires valid Fried, wavelength, and wind values" ); + } + tau_0( t0, m_lam_0 ); + } + + if( config.isSet( "atm.z_mean" ) ) + { + if( !math::isFinite( zm ) || zm < 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "configured mean layer height must be finite and nonnegative" ); + } + + const realT currentHeight = z_mean(); + if( currentHeight == 0 ) + { + if( zm != 0 ) + { + return internal::mxlib_error_report( + error_t::invalidconfig, + "a zero-height atmosphere cannot be rescaled to a positive mean height" ); + } + } + else + { + z_mean( zm ); + } + } + + return validate(); } extern template class aoAtmosphere; diff --git a/include/ao/analysis/aoSystem.hpp b/include/ao/analysis/aoSystem.hpp index 714c28a46..22d122140 100644 --- a/include/ao/analysis/aoSystem.hpp +++ b/include/ao/analysis/aoSystem.hpp @@ -1159,17 +1159,14 @@ class aoSystem */ iosT &dumpAOSystem( iosT &ios /**< [in] a std::ostream-like stream. */ ); - /// Setup the configurator to configure this class - /** - * todo: "\test Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]" - */ + /// Setup the configurator to configure this class. void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); - /// Load the configuration of this class from a configurator + /// Load the configuration of this class from a configurator. /** - * \todo: "\test Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]"" + * \returns `error_t::noerror` when the configured atmosphere is valid, or its typed validation error. */ - void loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); + error_t loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ ); }; template @@ -3056,7 +3053,7 @@ void aoSystem::setupConfig( app::appConfigurator &conf } template -void aoSystem::loadConfig( app::appConfigurator &config ) +error_t aoSystem::loadConfig( app::appConfigurator &config ) { // WFS if( config.isSet( "aosys.wfs" ) ) @@ -3238,8 +3235,10 @@ void aoSystem::loadConfig( app::appConfigurator &confi if( config.isSet( "aosys.starMag" ) ) starMag( smag ); - atm.loadConfig( config ); + const error_t atmosphereStatus = atm.loadConfig( config ); psd.loadConfig( config ); + + return atmosphereStatus; } extern template class aoSystem, std::ostream>; diff --git a/include/ao/analysis/clAOLinearPredictor.hpp b/include/ao/analysis/clAOLinearPredictor.hpp index 852bcffba..4f1e421fc 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. }; - std::vector m_PSDtn; ///< Working memory for the regularized PSD + /// 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. + }; - std::vector m_psd2s; ///< Working memory for the 2-sided regularized PSD + /// 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_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; -public: + std::vector m_regResults; ///< Per-scale telemetry collected when requested. - clAOLinearPredictor() - {} + regularizationReport m_regularizationReport; ///< Diagnostic summary of the latest search. + + 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; } @@ -182,19 +218,28 @@ struct clAOLinearPredictor go_lp.b( m_lp.m_c ); CLAOLP_BREADCRUMB; - realT ll = 0, ul = 0; - gmax_lp = go_lp.maxStableGain( ll, ul ); + rv = go_lp.maxStableGain( gmax_lp ); + if( rv != error_t::noerror ) + { + m_regularizationReport.status = regularizationStatus::calculationFailure; + return rv; + } if( gmax_lp > m_gmax_lp ) { gmax_lp = m_gmax_lp; } CLAOLP_BREADCRUMB; - gopt_lp = go_lp.optGainOpenLoop( var_lp, PSDt, PSDn, gmax_lp, false ); + rv = go_lp.optGainOpenLoop( gopt_lp, var_lp, PSDt, PSDn, gmax_lp, false ); + if( rv != error_t::noerror ) + { + m_regularizationReport.status = regularizationStatus::calculationFailure; + return rv; + } 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 +252,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 +267,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 +313,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 +325,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 +341,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; } @@ -295,30 +370,42 @@ struct clAOLinearPredictor go_lp.b( m_lp.m_c ); CLAOLP_BREADCRUMB; - realT ll = 0, ul = 0; - gmax_lp = go_lp.maxStableGain( ll, ul ); - gopt_lp = go_lp.optGainOpenLoop( var_lp, PSDt, PSDn, gmax_lp, false ); + rv = go_lp.maxStableGain( gmax_lp ); + if( rv != error_t::noerror ) + { + m_regularizationReport.status = regularizationStatus::calculationFailure; + return rv; + } + + rv = go_lp.optGainOpenLoop( gopt_lp, var_lp, PSDt, PSDn, gmax_lp, false ); + if( rv != error_t::noerror ) + { + m_regularizationReport.status = regularizationStatus::calculationFailure; + return rv; + } 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 +414,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 +432,7 @@ struct clAOLinearPredictor } } - return 0; + return error_t::noerror; } }; diff --git a/include/ao/analysis/clGainOpt.hpp b/include/ao/analysis/clGainOpt.hpp index 6f66d40e0..47cdb388f 100644 --- a/include/ao/analysis/clGainOpt.hpp +++ b/include/ao/analysis/clGainOpt.hpp @@ -28,9 +28,12 @@ #define clGainOpt_hpp #ifdef MX_INCLUDE_BOOST - #include +#include #endif +#include +#include +#include #include #include @@ -38,6 +41,8 @@ #include "../../sys/timeUtils.hpp" #include "../../math/constants.hpp" +#include "../../math/floatUtils.hpp" +#include "../../error/error_t.hpp" // #define ALLOW_F_ZERO @@ -65,6 +70,56 @@ struct clGainOpt typedef _realT realT; ///< The real data type typedef std::complex<_realT> complexT; ///< The complex data type + /// Termination state of a maximum-stable-gain search. + enum class maxStableGainStatus + { + notRun, ///< No search has been attempted. + crossingFound, ///< A qualifying Nyquist crossing was found. + invalidInput, ///< The frequency grid or derived Nyquist values were invalid. + noCrossing ///< No qualifying Nyquist crossing was found. + }; + + /// Diagnostic summary of a maximum-stable-gain search. + struct maxStableGainReport + { + maxStableGainStatus status{ maxStableGainStatus::notRun }; ///< Search termination state. + size_t lowerIndex{ std::numeric_limits::max() }; ///< Index below the selected crossing. + size_t upperIndex{ std::numeric_limits::max() }; ///< Index above the selected crossing. + realT lowerFrequency{ std::numeric_limits::quiet_NaN() }; ///< Frequency below the crossing. + realT upperFrequency{ std::numeric_limits::quiet_NaN() }; ///< Frequency above the crossing. + realT crossingFrequency{ std::numeric_limits::quiet_NaN() }; ///< Interpolated crossing frequency. + realT crossingReal{ std::numeric_limits::quiet_NaN() }; ///< Interpolated real Nyquist value. + realT gain{ std::numeric_limits::quiet_NaN() }; ///< Maximum stable gain at the crossing. + }; + + /// Termination state of an open-loop optimum-gain search. + enum class optGainStatus + { + notRun, ///< No search has been attempted. + converged, ///< The minimizer converged inside the search interval. + boundaryLimited, ///< The reported minimum lies on a search boundary. + invalidInput, ///< The PSDs, search controls, or requested interval were invalid. + stabilityFailure, ///< The automatic maximum-stable-gain search failed. + iterationLimit, ///< The minimizer exhausted its iteration limit. + calculationFailure ///< The minimizer threw or returned invalid output. + }; + + /// Diagnostic summary of an open-loop optimum-gain search. + struct optGainReport + { + optGainStatus status{ optGainStatus::notRun }; ///< Search termination state. + uintmax_t iterations{ 0 }; ///< Minimizer iterations attempted. + size_t evaluations{ 0 }; ///< Objective evaluations performed. + realT requestedMaximumGain{ std::numeric_limits::quiet_NaN() }; ///< Caller-supplied gain limit. + realT searchMinimumGain{ std::numeric_limits::quiet_NaN() }; ///< Final minimizer lower bound. + realT searchMaximumGain{ std::numeric_limits::quiet_NaN() }; ///< Final minimizer upper bound. + realT minimumEvaluatedGain{ std::numeric_limits::quiet_NaN() }; ///< Smallest evaluated gain. + realT maximumEvaluatedGain{ std::numeric_limits::quiet_NaN() }; ///< Largest evaluated gain. + realT gain{ std::numeric_limits::quiet_NaN() }; ///< Best gain returned by the minimizer. + realT variance{ std::numeric_limits::quiet_NaN() }; ///< Variance at the best gain. + maxStableGainReport stability; ///< Automatic stability-search diagnostics, when requested. + }; + protected: int m_N; ///< Number of integrations in the (optional) moving average. Default is 1. realT m_Ti; ///< The loop sampling interval @@ -76,9 +131,10 @@ struct clGainOpt std::vector m_f; ///< Vector of frequencies - bool m_fChanged{ true }; ///< True if frequency or max size of m_a and m_b changes + /// True when frequency, sampling interval, or required controller tap count invalidates m_cs and m_ss. + bool m_trigCacheChanged{ true }; - bool m_changed{ true }; ///< True if any of the members which make up the basic transfer functions are changed + bool m_changed{ true }; ///< True if any of the members which make up the basic transfer functions are changed Eigen::Array m_cs; Eigen::Array m_ss; @@ -170,8 +226,8 @@ struct clGainOpt /// Get a single FIR coefficient /** - * \returns a single FIR coefficient - */ + * \returns a single FIR coefficient + */ realT b( size_t i /**< [in] the index of the FIR coefficient*/ ) { return m_b[i]; @@ -199,8 +255,8 @@ struct clGainOpt coefficients, which is copied to m_a.*/ ); /// Get a single IIR coefficient /** - * \returns a single IIR coefficient - */ + * \returns a single IIR coefficient + */ realT a( size_t i ) { return m_a[i]; @@ -359,69 +415,40 @@ struct clGainOpt realT g ///< [in] the gain. ); - /// Find the maximum stable gain for the loop parameters - /** - * - * \returns the maximum stable gain for the loop parameters - */ - - /// Find the maximum stable gain for the loop parameters - /** Conducts a search along the Nyquist contour of the open-loop transfer function to find - * the most-negative crossing of the real axis. - * - * \returns the maximum stable gain for the loop parameters - */ - realT maxStableGain( realT &ll, ///< [in.out] the lower limit used for the search - realT &ul ///< [in.out] the upper limit used for hte search - ); - - /// Find the maximum stable gain for the loop parameters - /** Conducts a search along the Nyquist contour of the open-loop transfer function to find - * the most-negative crossing of the real axis. - * - * This version allows constant arguments. - * \overload - * - * \returns the maximum stable gain for the loop parameters - */ - realT maxStableGain( const realT &ll, ///< [in] the lower limit used for the search - const realT &ul ///< [in] the upper limit used for hte search - ); - /// Find the maximum stable gain for the loop parameters /** Conducts a search along the Nyquist contour of the open-loop transfer function to find * the most-negative crossing of the real axis. * - * This version uses m_maxFindMin for the lower limit and no upper limit. + * Crossings below m_maxFindMin are ignored. * - * \overload - * - * \returns the maximum stable gain for the loop parameters + * \returns `error_t::noerror` when a crossing is found, `error_t::notfound` when none is found, or an input error. */ - realT maxStableGain(); + mx::error_t maxStableGain( realT &gain, /**< [out] maximum stable gain; NaN on failure */ + maxStableGainReport *report = nullptr /**< [out] optional search diagnostics */ ); /// Return the optimum closed loop gain given an open loop PSD - /** Uses _gmax for the upper limit. - * \returns the optimum gain + /** Determines the maximum stable gain before minimizing the variance. + * \returns `error_t::noerror` on convergence or a boundary-limited result, otherwise an explicit failure status. */ - realT optGainOpenLoop( realT &var, ///< [out] the variance at the optimum gain - const std::vector &PSDerr, ///< [in] open loop error PSD - const std::vector &PSDnoise, ///< [in] open loop measurement noise PSD - bool gridSearch /**< [in] flag controlling whether an initial grid - search is done to find the global minimum*/ + mx::error_t optGainOpenLoop( realT &gain, ///< [out] optimum gain; NaN on failure + realT &var, ///< [out] variance at the optimum gain; NaN on failure + const std::vector &PSDerr, ///< [in] open-loop error PSD + const std::vector &PSDnoise, ///< [in] open-loop measurement-noise PSD + bool gridSearch, ///< [in] whether to perform a coarse initial search + optGainReport *report = nullptr ///< [out] optional search diagnostics ); /// Return the optimum closed loop gain given an open loop PSD /** - * \returns the optimum gain. + * \returns `error_t::noerror` on convergence or a boundary-limited result, otherwise an explicit failure status. */ - realT optGainOpenLoop( realT &var, ///< [out] the variance at the optimum gain - const std::vector &PSDerr, ///< [in] open loop error PSD - const std::vector &PSDnoise, ///< [in] open loop measurement noise PSD - realT &gmax, /**< [in] maximum gain to consider. - If 0, then _gmax is used.*/ - bool gridSearch /**< [in] flag controlling whether an initial grid - search is done to find the global minimum*/ + mx::error_t optGainOpenLoop( realT &gain, ///< [out] optimum gain; best estimate on timeout + realT &var, ///< [out] variance at the optimum gain + const std::vector &PSDerr, ///< [in] open-loop error PSD + const std::vector &PSDnoise, ///< [in] open-loop measurement-noise PSD + realT maximumGain, ///< [in] maximum stable gain bounding the search + bool gridSearch, ///< [in] whether to perform a coarse initial search + optGainReport *report = nullptr ///< [out] optional search diagnostics ); /// Calculate the pseudo open-loop PSD given a closed loop PSD @@ -468,7 +495,7 @@ void clGainOpt::init() m_minFindBits = std::numeric_limits::digits; m_minFindMaxIter = 10000; - m_fChanged = true; + m_trigCacheChanged = true; m_changed = true; } @@ -505,6 +532,7 @@ void clGainOpt::Ti( realT newTi ) } m_Ti = newTi; + m_trigCacheChanged = true; m_changed = true; } @@ -531,7 +559,7 @@ void clGainOpt::b( const std::vector &newB ) { if( newB.size() > (size_t)m_cs.cols() ) { - m_fChanged = true; + m_trigCacheChanged = true; } m_b = newB; @@ -543,7 +571,7 @@ void clGainOpt::b( const Eigen::Array &newB ) { if( newB.cols() > m_cs.cols() ) { - m_fChanged = true; + m_trigCacheChanged = true; } m_b.resize( newB.cols() ); @@ -558,21 +586,21 @@ void clGainOpt::b( const Eigen::Array &newB ) template void clGainOpt::bScale( realT scale ) +{ + for( size_t n = 0; n < m_b.size(); ++n ) { - for( size_t n = 0; n < m_b.size(); ++n ) - { - m_b[n] *= scale; - } - - m_changed = true; + m_b[n] *= scale; } + m_changed = true; +} + template void clGainOpt::a( const std::vector &newA ) { if( newA.size() + 1 > (size_t)m_cs.cols() ) { - m_fChanged = true; + m_trigCacheChanged = true; } m_a = newA; @@ -584,7 +612,7 @@ void clGainOpt::a( const Eigen::Array &newA ) { if( newA.cols() + 1 > m_cs.cols() ) { - m_fChanged = true; + m_trigCacheChanged = true; } m_a.resize( newA.cols() ); @@ -611,14 +639,12 @@ void clGainOpt::aScale( realT scale ) template void clGainOpt::remember( const realT &rem ) { - if(m_remember != rem) + if( m_remember != rem ) { m_remember = rem; m_changed = true; } - - } template @@ -630,20 +656,20 @@ realT clGainOpt::remember() template void clGainOpt::setLeakyIntegrator( realT remember ) { - if(m_b.size() != 1 || m_a.size() != 1 || m_b[0] != 1.0 || m_a[0] != 1.0 || m_remember != remember) + if( m_b.size() != 1 || m_a.size() != 1 || m_b[0] != 1.0 || m_a[0] != 1.0 || m_remember != remember ) { - if(m_b.size() != 1) + if( m_b.size() != 1 ) { m_b.resize( 1 ); - m_fChanged = true; + m_trigCacheChanged = true; } m_b[0] = 1.0; - if(m_a.size() != 1) + if( m_a.size() != 1 ) { m_a.resize( 1 ); - m_fChanged = true; + m_trigCacheChanged = true; } m_a[0] = 1.0; @@ -663,7 +689,7 @@ void clGainOpt::f( realT *newF, size_t nF ) m_f[i] = newF[i]; } - m_fChanged = true; + m_trigCacheChanged = true; m_changed = true; } @@ -671,7 +697,7 @@ template void clGainOpt::f( const std::vector &newF ) { m_f = newF; - m_fChanged = true; + m_trigCacheChanged = true; m_changed = true; } @@ -716,7 +742,7 @@ std::complex clGainOpt::olXfer( int fi, complexT &H_dm, complexT & } #ifdef PRECALC_TRIG - if( m_fChanged ) + if( m_trigCacheChanged ) { size_t jmax = std::max( m_a.size() + 1, m_b.size() ); @@ -735,7 +761,7 @@ std::complex clGainOpt::olXfer( int fi, complexT &H_dm, complexT & } } - m_fChanged = false; + m_trigCacheChanged = false; } #endif @@ -1003,51 +1029,81 @@ realT clGainOpt::clVariance( const std::vector &PSDerr, const std: } template -realT clGainOpt::maxStableGain( realT &ll, realT &ul ) +mx::error_t clGainOpt::maxStableGain( realT &gain, maxStableGainReport *report ) { - static_cast( ul ); + maxStableGainReport localReport; + maxStableGainReport &activeReport = report == nullptr ? localReport : *report; + activeReport = {}; + gain = std::numeric_limits::quiet_NaN(); - std::vector re, im; + if( m_f.size() < 2 ) + { + activeReport.status = maxStableGainStatus::invalidInput; + return error_t::sizeerr; + } - if( ll == 0 ) - ll = m_maxFindMin; + for( size_t index = 0; index < m_f.size(); ++index ) + { + if( !math::isFinite( m_f[index] ) || m_f[index] < 0 || ( index > 0 && m_f[index] <= m_f[index - 1] ) ) + { + activeReport.status = maxStableGainStatus::invalidInput; + return error_t::invalidarg; + } + } + + std::vector re, im; nyquist( re, im, 1.0 ); - int gi_c = re.size() - 1; + for( size_t index = 0; index < re.size(); ++index ) + { + if( !math::isFinite( re[index] ) || !math::isFinite( im[index] ) ) + { + activeReport.status = maxStableGainStatus::invalidInput; + return error_t::error; + } + } - for( int gi = re.size() - 2; gi >= 0; --gi ) + bool crossingFound = false; + for( size_t index = 0; index + 1 < re.size(); ++index ) { - if( -1.0 / re[gi] < ll ) + if( !( im[index] < 0 && im[index + 1] >= 0 ) ) + { continue; + } - if( ( re[gi] < 0 ) && ( im[gi + 1] >= 0 && im[gi] < 0 ) ) + const realT fraction = -im[index] / ( im[index + 1] - im[index] ); + const realT crossingReal = re[index] + fraction * ( re[index + 1] - re[index] ); + const realT crossingGain = -realT( 1 ) / crossingReal; + if( crossingReal >= 0 || !math::isFinite( crossingGain ) || crossingGain < m_maxFindMin ) { - // Check for loop back in Nyquist diagram - if( re[gi] <= re[gi_c] ) - gi_c = gi; + continue; } - } - - return -1.0 / re[gi_c]; -} -template -realT maxStableGain( const realT &ll, const realT &ul ) -{ - realT rll = ll; - realT rul = ul; + if( crossingFound && crossingReal >= activeReport.crossingReal ) + { + continue; + } - maxStableGain( rll, rul ); -} + crossingFound = true; + activeReport.lowerIndex = index; + activeReport.upperIndex = index + 1; + activeReport.lowerFrequency = m_f[index]; + activeReport.upperFrequency = m_f[index + 1]; + activeReport.crossingFrequency = m_f[index] + fraction * ( m_f[index + 1] - m_f[index] ); + activeReport.crossingReal = crossingReal; + activeReport.gain = crossingGain; + } -template -realT clGainOpt::maxStableGain() -{ - realT ul = 0; - realT ll = m_maxFindMin; + if( !crossingFound ) + { + activeReport.status = maxStableGainStatus::noCrossing; + return error_t::notfound; + } - return maxStableGain( ll, ul ); + activeReport.status = maxStableGainStatus::crossingFound; + gain = activeReport.gain; + return error_t::noerror; } // Implement the minimization, allowing pre-compiled specializations @@ -1055,158 +1111,263 @@ namespace impl { template -realT optGainOpenLoop( clGainOptOptGain_OL &olgo, - realT &var, - const realT &gmax, - const realT &minFindMin, - const realT &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ) +/// Minimize an open-loop variance objective on a bounded gain interval. +mx::error_t optGainOpenLoop( realT &gain, ///< [out] best gain estimate + realT &var, ///< [out] variance at the best gain estimate + clGainOptOptGain_OL &olgo, ///< [in,out] variance objective and diagnostics + const realT &minimumGain, ///< [in] lower gain bound + const realT &maximumGain, ///< [in] upper gain bound + int minFindBits, ///< [in] requested precision in binary digits + uintmax_t minFindMaxIter, ///< [in] maximum minimizer iterations + uintmax_t &iters ///< [out] minimizer iterations used +) { #ifdef MX_INCLUDE_BOOST - realT gopt; + gain = std::numeric_limits::quiet_NaN(); + var = std::numeric_limits::quiet_NaN(); try { std::pair brack; brack = boost::math::tools::brentm_findm_minima, realT>( olgo, - minFindMin, - minFindMaxFact * gmax, + minimumGain, + maximumGain, minFindBits, minFindMaxIter, iters ); - gopt = brack.first; + gain = brack.first; var = brack.second; } catch( ... ) { - std::cerr << "optGainOpenLoop: No root found\n"; - gopt = minFindMaxFact * gmax; - var = 0; + return error_t::exception; + } + + if( iters >= minFindMaxIter ) + { + return error_t::timeout; } - return gopt; + return error_t::noerror; #else static_assert( std::is_fundamental::value || !std::is_fundamental::value, - "impl::optGainOpenLoop is not specialized for type realT, and MX_INCLUDE_BOOST is not " - "defined, so I can't just use boost." ); - return 0; + "impl::optGainOpenLoop is not specialized for type realT, and MX_INCLUDE_BOOST is not " + "defined, so I can't just use boost." ); + return error_t::notimpl; #endif } template <> -float optGainOpenLoop( clGainOptOptGain_OL &olgo, - float &var, - const float &gmax, - const float &minFindMin, - const float &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ); +/// Float specialization of the bounded open-loop gain minimizer. +mx::error_t optGainOpenLoop( float &gain, ///< [out] best gain estimate + float &var, ///< [out] variance at the best gain estimate + clGainOptOptGain_OL &olgo, ///< [in,out] variance objective and diagnostics + const float &minimumGain, ///< [in] lower gain bound + const float &maximumGain, ///< [in] upper gain bound + int minFindBits, ///< [in] requested precision in binary digits + uintmax_t minFindMaxIter, ///< [in] maximum minimizer iterations + uintmax_t &iters ///< [out] minimizer iterations used +); template <> -double optGainOpenLoop( clGainOptOptGain_OL &olgo, - double &var, - const double &gmax, - const double &minFindMin, - const double &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ); +/// Double specialization of the bounded open-loop gain minimizer. +mx::error_t optGainOpenLoop( double &gain, ///< [out] best gain estimate + double &var, ///< [out] variance at the best gain estimate + clGainOptOptGain_OL &olgo, ///< [in,out] variance objective and diagnostics + const double &minimumGain, ///< [in] lower gain bound + const double &maximumGain, ///< [in] upper gain bound + int minFindBits, ///< [in] requested precision in binary digits + uintmax_t minFindMaxIter, ///< [in] maximum minimizer iterations + uintmax_t &iters ///< [out] minimizer iterations used +); template <> -long double optGainOpenLoop( clGainOptOptGain_OL &olgo, - long double &var, - const long double &gmax, - const long double &minFindMin, - const long double &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ); +/// Long-double specialization of the bounded open-loop gain minimizer. +mx::error_t +optGainOpenLoop( long double &gain, ///< [out] best gain estimate + long double &var, ///< [out] variance at the best gain estimate + clGainOptOptGain_OL &olgo, ///< [in,out] variance objective and diagnostics + const long double &minimumGain, ///< [in] lower gain bound + const long double &maximumGain, ///< [in] upper gain bound + int minFindBits, ///< [in] requested precision in binary digits + uintmax_t minFindMaxIter, ///< [in] maximum minimizer iterations + uintmax_t &iters ///< [out] minimizer iterations used +); #ifdef HASQUAD template <> -_m_float128 optGainOpenLoop<_m_float128>( clGainOptOptGain_OL<_m_float128> &olgo, - _m_float128 &var, - const _m_float128 &gmax, - const _m_float128 &minFindMin, - const _m_float128 &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ); +/// Quad-precision specialization of the bounded open-loop gain minimizer. +mx::error_t +optGainOpenLoop<_m_float128>( _m_float128 &gain, ///< [out] best gain estimate + _m_float128 &var, ///< [out] variance at the best gain estimate + clGainOptOptGain_OL<_m_float128> &olgo, ///< [in,out] variance objective and diagnostics + const _m_float128 &minimumGain, ///< [in] lower gain bound + const _m_float128 &maximumGain, ///< [in] upper gain bound + int minFindBits, ///< [in] requested precision in binary digits + uintmax_t minFindMaxIter, ///< [in] maximum minimizer iterations + uintmax_t &iters ///< [out] minimizer iterations used +); #endif } // namespace impl template -realT clGainOpt::optGainOpenLoop( realT &var, - const std::vector &PSDerr, - const std::vector &PSDnoise, - bool gridSearch ) +mx::error_t clGainOpt::optGainOpenLoop( realT &gain, + realT &var, + const std::vector &PSDerr, + const std::vector &PSDnoise, + bool gridSearch, + optGainReport *report ) { - realT gmax = 0; - return optGainOpenLoop( var, PSDerr, PSDnoise, gmax, gridSearch ); + maxStableGainReport stabilityReport; + realT maximumGain; + error_t rv = maxStableGain( maximumGain, &stabilityReport ); + if( rv != error_t::noerror ) + { + gain = std::numeric_limits::quiet_NaN(); + var = std::numeric_limits::quiet_NaN(); + if( report != nullptr ) + { + *report = {}; + report->status = optGainStatus::stabilityFailure; + report->stability = stabilityReport; + } + return rv; + } + + optGainReport optimizationReport; + rv = optGainOpenLoop( gain, var, PSDerr, PSDnoise, maximumGain, gridSearch, &optimizationReport ); + optimizationReport.stability = stabilityReport; + if( report != nullptr ) + { + *report = optimizationReport; + } + return rv; } template -realT clGainOpt::optGainOpenLoop( - realT &var, const std::vector &PSDerr, const std::vector &PSDnoise, realT &gmax, bool gridSearch ) +mx::error_t clGainOpt::optGainOpenLoop( realT &gain, + realT &var, + const std::vector &PSDerr, + const std::vector &PSDnoise, + realT maximumGain, + bool gridSearch, + optGainReport *report ) { + optGainReport localReport; + optGainReport &activeReport = report == nullptr ? localReport : *report; + activeReport = {}; + activeReport.requestedMaximumGain = maximumGain; + gain = std::numeric_limits::quiet_NaN(); + var = std::numeric_limits::quiet_NaN(); + + if( m_f.size() < 2 || PSDerr.size() != m_f.size() || PSDnoise.size() != m_f.size() ) + { + activeReport.status = optGainStatus::invalidInput; + return error_t::sizeerr; + } + + if( !math::isFinite( maximumGain ) || maximumGain <= 0 || !math::isFinite( m_minFindMin ) || m_minFindMin < 0 || + !math::isFinite( m_minFindMaxFact ) || m_minFindMaxFact <= 0 || m_minFindMaxFact > 1 || m_minFindBits <= 0 || + m_minFindMaxIter == 0 ) + { + activeReport.status = optGainStatus::invalidInput; + return error_t::invalidconfig; + } + + const realT requestedMinimum = m_minFindMin; + const realT requestedMaximum = m_minFindMaxFact * maximumGain; + if( !math::isFinite( requestedMaximum ) || requestedMaximum <= requestedMinimum ) + { + activeReport.status = optGainStatus::invalidInput; + return error_t::invalidconfig; + } + clGainOptOptGain_OL olgo; olgo.go = this; olgo.PSDerr = &PSDerr; olgo.PSDnoise = &PSDnoise; - if( gmax <= 0 ) - { - gmax = maxStableGain(); - } - - realT ming = m_minFindMin; - realT maxg = gmax; + realT minimumGain = requestedMinimum; + realT maximumSearchGain = requestedMaximum; + bool searchBoundarySelected = false; if( gridSearch ) { - realT gstpsz = 0.05; - realT gg = m_minFindMaxFact * gmax; - realT var0 = clVariance( PSDerr, PSDnoise, gg ); - realT mingg = gg; + const realT gainStep = std::min( realT( 0.05 ), requestedMaximum - requestedMinimum ); + realT currentGain = requestedMaximum; + realT minimumVariance = olgo( currentGain ); + realT gainAtMinimum = currentGain; - while( gg > m_minFindMin ) + while( currentGain > requestedMinimum ) { - gg -= gstpsz; - realT var1 = clVariance( PSDerr, PSDnoise, gg ); + const realT nextGain = std::max( requestedMinimum, currentGain - gainStep ); + if( nextGain >= currentGain ) + { + break; + } + + currentGain = nextGain; + const realT candidateVariance = olgo( currentGain ); - if( var1 < var0 ) + if( candidateVariance < minimumVariance ) { - var0 = var1; - mingg = gg; + minimumVariance = candidateVariance; + gainAtMinimum = currentGain; } } - ming = mingg - gstpsz; - maxg = mingg + gstpsz; + minimumGain = std::max( requestedMinimum, gainAtMinimum - gainStep ); + maximumSearchGain = std::min( requestedMaximum, gainAtMinimum + gainStep ); + searchBoundarySelected = gainAtMinimum == requestedMinimum || gainAtMinimum == requestedMaximum; + } - if( ming < m_minFindMin ) - ming = m_minFindMin; - if( maxg > gmax ) - maxg = gmax; + activeReport.searchMinimumGain = minimumGain; + activeReport.searchMaximumGain = maximumSearchGain; + + uintmax_t iterations = m_minFindMaxIter; + error_t rv = impl::optGainOpenLoop( gain, + var, + olgo, + minimumGain, + maximumSearchGain, + m_minFindBits, + m_minFindMaxIter, + iterations ); + + activeReport.iterations = iterations; + activeReport.evaluations = olgo.evaluations; + activeReport.minimumEvaluatedGain = olgo.minimumEvaluatedGain; + activeReport.maximumEvaluatedGain = olgo.maximumEvaluatedGain; + activeReport.gain = gain; + activeReport.variance = var; + + if( rv == error_t::timeout ) + { + activeReport.status = optGainStatus::iterationLimit; + return rv; } - uintmax_t iters; - realT val = - impl::optGainOpenLoop( olgo, var, maxg, ming, m_minFindMaxFact, m_minFindBits, m_minFindMaxIter, iters ); + if( rv != error_t::noerror || !math::isFinite( gain ) || !math::isFinite( var ) ) + { + activeReport.status = optGainStatus::calculationFailure; + gain = std::numeric_limits::quiet_NaN(); + var = std::numeric_limits::quiet_NaN(); + activeReport.gain = gain; + activeReport.variance = var; + return rv == error_t::noerror ? error_t::error : rv; + } - if( iters >= m_minFindMaxIter ) + if( searchBoundarySelected || gain <= minimumGain || gain >= maximumSearchGain ) { - // #pragma omp critical - { - std::cerr << "\nclGainOpt::optGainOpenLoop: minFindMaxIter (" << m_minFindMaxIter << ") reached\n"; - } + activeReport.status = optGainStatus::boundaryLimited; + } + else + { + activeReport.status = optGainStatus::converged; } - return val; + return error_t::noerror; } template @@ -1248,12 +1409,19 @@ int clGainOpt::nyquist( std::vector &re, std::vector &im, r template struct clGainOptOptGain_OL { - clGainOpt *go; - const std::vector *PSDerr; - const std::vector *PSDnoise; - - realT operator()( const realT &g ) + clGainOpt *go{ nullptr }; ///< Gain optimizer used to evaluate variance. + const std::vector *PSDerr{ nullptr }; ///< Open-loop disturbance PSD. + const std::vector *PSDnoise{ nullptr }; ///< Measurement-noise PSD. + size_t evaluations{ 0 }; ///< Objective evaluations performed. + realT minimumEvaluatedGain{ std::numeric_limits::max() }; ///< Smallest gain evaluated. + realT maximumEvaluatedGain{ std::numeric_limits::lowest() }; ///< Largest gain evaluated. + + /// Evaluate closed-loop variance at a candidate gain and update diagnostics. + realT operator()( const realT &g /**< [in] candidate gain */ ) { + ++evaluations; + minimumEvaluatedGain = std::min( minimumEvaluatedGain, g ); + maximumEvaluatedGain = std::max( maximumEvaluatedGain, g ); return go->clVariance( *PSDerr, *PSDnoise, g ); } }; diff --git a/include/ao/analysis/fourierTemporalPSD.hpp b/include/ao/analysis/fourierTemporalPSD.hpp index de0aea1fd..f9d3a3914 100644 --- a/include/ao/analysis/fourierTemporalPSD.hpp +++ b/include/ao/analysis/fourierTemporalPSD.hpp @@ -27,8 +27,15 @@ #ifndef fourierTemporalPSD_hpp #define fourierTemporalPSD_hpp +#include +#include +#include #include #include +#include +#include +#include +#include #include @@ -39,6 +46,7 @@ #include "../../mxlib.hpp" #include "../../math/constants.hpp" +#include "../../math/floatUtils.hpp" #include "../../math/func/jinc.hpp" #include "../../math/func/airyPattern.hpp" #include "../../math/vectorUtils.hpp" @@ -71,19 +79,260 @@ namespace analysis #ifndef WSZ - /** \def WFZ - * \brief Size of the GSL integration workspace - */ - #define WSZ 100000 +/** \def WFZ + * \brief Size of the GSL integration workspace + */ +#define WSZ 100000 #endif enum basis : unsigned int { - basic, ///< The basic sine and cosine Fourier modes + basic, ///< The basic sine and cosine Fourier modes modified ///< The modified Fourier basis from \cite males_guyon_2017 }; +/// Policy for handling GSL quadrature non-convergence statuses. +enum class fourierTemporalPSDPolicy +{ + permissive, ///< Retain the best finite approximation and record the status. + strict ///< Record every status and return an error if any integration does not converge. +}; + +/// Aggregated GSL quadrature diagnostics for a Fourier temporal PSD calculation. +template +struct fourierTemporalPSDReport +{ + /// Summary of one GSL status code. + struct statusSummary + { + size_t count{ 0 }; ///< Number of occurrences of this status. + std::map countByLayer; ///< Number of occurrences in each atmospheric layer. + realT maximumAbsoluteError{ 0 }; ///< Largest GSL absolute-error estimate. + realT maximumToleranceRatio{ 0 }; ///< Largest error estimate relative to the requested tolerance. + size_t worstLayer{ 0 }; ///< Layer containing the largest tolerance ratio. + realT worstFrequency{ 0 }; ///< Frequency containing the largest tolerance ratio. + }; + + size_t integrationsAttempted{ 0 }; ///< Total number of quadrature calls. + size_t integrationsConverged{ 0 }; ///< Number of quadrature calls returning `GSL_SUCCESS`. + std::map gslStatus; ///< Summaries keyed by the raw GSL status code. + + /// Reset all accumulated diagnostics. + void clear(); + + /// Record one quadrature result. + void record( int status, /**< [in] raw GSL status code */ + size_t layer, /**< [in] atmospheric layer index */ + realT frequency, /**< [in] temporal frequency */ + realT result, /**< [in] quadrature result */ + realT absoluteError, /**< [in] GSL absolute-error estimate */ + realT absoluteTolerance, /**< [in] requested absolute tolerance */ + realT relativeTolerance /**< [in] requested relative tolerance */ ); + + /// Merge another report into this report. + void merge( const fourierTemporalPSDReport &other /**< [in] report to merge */ ); + + /// Return the total number of non-successful integrations. + [[nodiscard]] size_t failureCount() const; + + /// Write a human-readable summary of the accumulated quadrature diagnostics. + void write( std::ostream &output /**< [out] stream receiving the summary */ ) const; +}; + +/// \cond fourierTemporalPSD_detail +namespace fourierTemporalPSD_detail +{ + +/// Function type used to allocate a GSL integration workspace. +using gslWorkspaceAllocator = gsl_integration_workspace *(*)( size_t ); + +/// Deleter providing RAII ownership for a GSL integration workspace. +struct gslWorkspaceDeleter +{ + /// Free an allocated GSL integration workspace. + void operator()( gsl_integration_workspace *workspace /**< [in] workspace to free */ ) const noexcept + { + if( workspace != nullptr ) + { + gsl_integration_workspace_free( workspace ); + } + } +}; + +/// Unique ownership handle for a GSL integration workspace. +using gslWorkspacePtr = std::unique_ptr; + +/// Return whether a GSL status represents a potentially usable non-converged approximation. +inline bool isConvergenceStatus( int status ) +{ + return status == GSL_EMAXITER || status == GSL_EROUND || status == GSL_ESING || status == GSL_EDIVERGE; +} + +/// Convert a fatal GSL status to an mxlib status. +inline error_t gslStatusToError( int status ) +{ + if( status == GSL_ENOMEM ) + { + return error_t::allocerr; + } + + if( status == GSL_EDOM || status == GSL_EINVAL ) + { + return error_t::invalidconfig; + } + + return error_t::liberr; +} + +/// Apply the configured non-convergence policy to one GSL status. +inline error_t applyPolicy( int status, fourierTemporalPSDPolicy policy ) +{ + if( status == GSL_SUCCESS ) + { + return error_t::noerror; + } + + if( isConvergenceStatus( status ) ) + { + return policy == fourierTemporalPSDPolicy::permissive ? error_t::noerror : error_t::liberr; + } + + return gslStatusToError( status ); +} + +/// Mutex serializing scoped changes to GSL's process-global error handler. +inline std::mutex &gslErrorHandlerMutex() +{ + static std::mutex mutex; + return mutex; +} + +/// Disable the GSL error handler for a complete top-level PSD calculation and restore it on exit. +/** The mutex serializes handler changes made by this implementation. Unrelated code cannot be protected from GSL's + * process-global handler state unless it coordinates with the same mutex. + */ +class scopedGslErrorHandlerOff +{ + public: + /// Lock handler management and retain the previously installed handler. + scopedGslErrorHandlerOff() : m_lock( gslErrorHandlerMutex() ), m_previous( gsl_set_error_handler_off() ) + { + } + + /// Disallow copying ownership of the saved handler. + scopedGslErrorHandlerOff( const scopedGslErrorHandlerOff & ) = delete; + + /// Disallow copy assignment of the handler guard. + scopedGslErrorHandlerOff &operator=( const scopedGslErrorHandlerOff & ) = delete; + + /// Restore the previously installed handler before releasing the lock. + ~scopedGslErrorHandlerOff() + { + static_cast( gsl_set_error_handler( m_previous ) ); + } + + private: + std::unique_lock m_lock; ///< Lock held while the handler is disabled. + gsl_error_handler_t *m_previous{ nullptr }; ///< Handler restored on destruction. +}; + +} // namespace fourierTemporalPSD_detail +/// \endcond + +template +void fourierTemporalPSDReport::clear() +{ + integrationsAttempted = 0; + integrationsConverged = 0; + gslStatus.clear(); +} + +template +void fourierTemporalPSDReport::record( int status, + size_t layer, + realT frequency, + realT result, + realT absoluteError, + realT absoluteTolerance, + realT relativeTolerance ) +{ + ++integrationsAttempted; + if( status == GSL_SUCCESS ) + { + ++integrationsConverged; + return; + } + + statusSummary &summary = gslStatus[status]; + ++summary.count; + ++summary.countByLayer[layer]; + summary.maximumAbsoluteError = std::max( summary.maximumAbsoluteError, std::abs( absoluteError ) ); + + const realT requestedTolerance = std::max( std::abs( absoluteTolerance ), std::abs( relativeTolerance * result ) ); + const realT toleranceRatio = requestedTolerance > 0 ? std::abs( absoluteError ) / requestedTolerance + : std::numeric_limits::infinity(); + if( toleranceRatio >= summary.maximumToleranceRatio ) + { + summary.maximumToleranceRatio = toleranceRatio; + summary.worstLayer = layer; + summary.worstFrequency = frequency; + } +} + +template +void fourierTemporalPSDReport::merge( const fourierTemporalPSDReport &other ) +{ + integrationsAttempted += other.integrationsAttempted; + integrationsConverged += other.integrationsConverged; + + for( const auto &[status, otherSummary] : other.gslStatus ) + { + statusSummary &summary = gslStatus[status]; + summary.count += otherSummary.count; + for( const auto &[layer, count] : otherSummary.countByLayer ) + { + summary.countByLayer[layer] += count; + } + summary.maximumAbsoluteError = std::max( summary.maximumAbsoluteError, otherSummary.maximumAbsoluteError ); + if( otherSummary.maximumToleranceRatio >= summary.maximumToleranceRatio ) + { + summary.maximumToleranceRatio = otherSummary.maximumToleranceRatio; + summary.worstLayer = otherSummary.worstLayer; + summary.worstFrequency = otherSummary.worstFrequency; + } + } +} + +template +size_t fourierTemporalPSDReport::failureCount() const +{ + return integrationsAttempted - integrationsConverged; +} + +template +void fourierTemporalPSDReport::write( std::ostream &output ) const +{ + output << "GSL quadrature: " << integrationsConverged << '/' << integrationsAttempted << " converged\n"; + for( const auto &[status, summary] : gslStatus ) + { + output << " " << gsl_strerror( status ) << " (" << status << "): " << summary.count << ", max absolute error " + << summary.maximumAbsoluteError << ", max tolerance ratio " << summary.maximumToleranceRatio + << " at layer " << summary.worstLayer << ", frequency " << summary.worstFrequency << ", layers {"; + bool firstLayer = true; + for( const auto &[layer, count] : summary.countByLayer ) + { + if( !firstLayer ) + { + output << ", "; + } + output << layer << ": " << count; + firstLayer = false; + } + output << "}\n"; + } +} + // Forward declaration template realT F_basic( realT kv, void *params ); @@ -92,7 +341,6 @@ realT F_basic( realT kv, void *params ); template realT F_mod( realT kv, void *params ); - /// Class to manage the calculation of temporal PSDs of the Fourier modes in atmospheric turbulence. /** Works with both basic (sines/cosines) and modified Fourier modes. * @@ -101,8 +349,6 @@ realT F_mod( realT kv, void *params ); * * \todo Split off the integration parameters in a separate structure. * \todo once integration parameters are in a separate structure, make this a class with protected members. - * \todo GSL error handling - * * \ingroup mxAOAnalytic */ template @@ -114,6 +360,9 @@ struct fourierTemporalPSD /// The complex type for arithmetic typedef std::complex complexT; + /// Quadrature report type used by this specialization. + typedef fourierTemporalPSDReport reportT; + /// Pointer to an AO system structure. aosysT *m_aosys{ nullptr }; @@ -136,9 +385,14 @@ struct fourierTemporalPSD int _useBasis; ///< Set to MXAO_FTPSD_BASIS_BASIC/MODIFIED/PROJECTED_* to use the basic sin/cos modes, the modified ///< Fourier modes, or a projection of them. - /// Workspace for the gsl integrators, allocated to WSZ if constructed as worker (with allocate == true). - gsl_integration_workspace *_w; + protected: + /// Unique ownership of the GSL integration workspace used by worker instances. + fourierTemporalPSD_detail::gslWorkspacePtr m_workspace; + + /// Allocation function used when a worker lazily creates its GSL workspace. + fourierTemporalPSD_detail::gslWorkspaceAllocator m_workspaceAllocator{ gsl_integration_workspace_alloc }; + public: realT _absTol; ///< The absolute tolerance to use in the GSL integrator realT _relTol; ///< The relative tolerance to use in the GSL integrator @@ -175,21 +429,45 @@ struct fourierTemporalPSD /// Default c'tor fourierTemporalPSD(); - /// Constructor with workspace allocation - /** - * \param allocate if true, then the workspace for GSL integrators is allocated. - */ - explicit fourierTemporalPSD( bool allocate ); + /// Disallow copying unique workspace ownership. + fourierTemporalPSD( const fourierTemporalPSD & ) = delete; - /// Destructor - /** Frees GSL workspace if it was allocated. - */ - ~fourierTemporalPSD(); + /// Disallow copy assignment of unique workspace ownership. + fourierTemporalPSD &operator=( const fourierTemporalPSD & ) = delete; + + /// Move workspace ownership and evaluator state. + fourierTemporalPSD( fourierTemporalPSD && ) noexcept = default; + + /// Move-assign workspace ownership and evaluator state. + fourierTemporalPSD &operator=( fourierTemporalPSD && ) noexcept = default; + + /// Release owned resources. + ~fourierTemporalPSD() = default; protected: + /// Construct with a custom workspace allocator. + explicit fourierTemporalPSD( + fourierTemporalPSD_detail::gslWorkspaceAllocator allocator /**< [in] workspace allocation function */ ); + /// Initialize parameters to default values. void initialize(); + /// Allocate the worker workspace if it is not already available. + error_t allocateWorkspace(); + + /// Validate state and arguments shared by single- and multilayer calculations. + error_t validatePsdInputs( const std::vector &PSD, /**< [in] output storage to validate */ + const std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax, /**< [in] maximum exactly integrated frequency */ + int layer_i, /**< [in] layer index, or -1 to validate all layers */ + fourierTemporalPSDPolicy policy /**< [in] non-convergence policy */ ); + + /// Validate the configured atmosphere and optionally a requested layer. + error_t validateAtmosphere( int layer_i /**< [in] layer index, or -1 to validate all layers */ ); + public: /** \name GSL Integration Tolerances * For good results it seems that absolute tolerance (absTol) needs to be 1e-10. Lower tolerances cause some @@ -236,25 +514,45 @@ struct fourierTemporalPSD */ realT fastestPeak( int m, int n ); + protected: + /// Calculate a single-layer temporal PSD while the caller manages the GSL error handler. + error_t singleLayerPSDImpl( std::vector &PSD, /**< [out] calculated PSD */ + std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int layer_i, /**< [in] atmospheric-layer index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax, /**< [in] maximum exactly integrated frequency */ + reportT &report, /**< [out] accumulated quadrature report */ + fourierTemporalPSDPolicy policy /**< [in] non-convergence policy */ ); + + public: /// Calculate the temporal PSD for a Fourier mode for a single layer. - /** + /** `PSD` and `freq` must have the same nonzero size. Frequencies must be finite, nonnegative, and strictly + * increasing. The AO system, integration controls, requested layer, and atmosphere are validated before + * calculation. A precondition or allocation failure leaves `PSD` unchanged and clears `report` when supplied. * - * \todo implement error checking. - * \todo need a way to track convergence failures in integral without throwing an error. - * \todo need better handling of averaging for the -17/3 extension. + * When extending beyond `fmax`, up to the last 50 exactly integrated bins are averaged after projection to the + * first tail frequency. If fewer than 50 exact bins are available, all available exact bins are used. At least one + * exact bin is required to initialize the tail. * + * In permissive mode, finite best approximations returned with `GSL_EMAXITER`, `GSL_EROUND`, `GSL_ESING`, or + * `GSL_EDIVERGE` are retained and summarized in `report`. In strict mode the calculation continues to characterize + * all such failures but returns `error_t::liberr` and the output PSD must be discarded. + * + * \returns `error_t::noerror` on success, an argument/configuration error for invalid inputs, or + * `error_t::liberr` when strict quadrature handling detects non-convergence. */ - int - singleLayerPSD( std::vector &PSD, ///< [out] the calculated PSD - std::vector &freq, ///< [in] the populated temporal frequency grid defining the frequencies - ///< at which the PSD is calculated - realT m, ///< [in] the first index of the spatial frequency - realT n, ///< [in] the second index of the spatial frequency - int layer_i, ///< [in] the index of the layer, for accessing the atmosphere parameters - int p, ///< [in] sets which mode is calculated (if basic modes, p = -1 for sine, p = +1 for cosine) - realT fmax = 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The PSD - ///< is filled in with a -17/3 power law past this frequency. - ); + error_t singleLayerPSD( + std::vector &PSD, /**< [out] calculated PSD */ + std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int layer_i, /**< [in] atmospheric-layer index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax = 0, /**< [in] maximum exactly integrated frequency, or 0 for the grid maximum */ + reportT *report = nullptr, /**< [out] optional quadrature report */ + fourierTemporalPSDPolicy policy = fourierTemporalPSDPolicy::permissive /**< [in] non-convergence policy */ ); ///\cond multilayerm_parallel // Conditional to exclude from Doxygen. @@ -267,48 +565,54 @@ struct fourierTemporalPSD }; // Parallelized version of multiLayerPSD, with OMP directives. - int m_multiLayerPSD( std::vector &PSD, - std::vector &freq, - realT m, - realT n, - int p, - realT fmax, - isParallel parallel ); + error_t multiLayerPSDImpl( std::vector &PSD, /**< [out] calculated PSD */ + std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax, /**< [in] maximum exactly integrated frequency */ + reportT &report, /**< [out] accumulated quadrature report */ + fourierTemporalPSDPolicy policy, /**< [in] non-convergence policy */ + isParallel parallel /**< [in] parallel dispatch tag */ ); // Non-Parallelized version of multiLayerPSD, without OMP directives. - int m_multiLayerPSD( std::vector &PSD, - std::vector &freq, - realT m, - realT n, - int p, - realT fmax, - isParallel parallel ); + error_t multiLayerPSDImpl( std::vector &PSD, /**< [out] calculated PSD */ + std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax, /**< [in] maximum exactly integrated frequency */ + reportT &report, /**< [out] accumulated quadrature report */ + fourierTemporalPSDPolicy policy, /**< [in] non-convergence policy */ + isParallel parallel /**< [in] sequential dispatch tag */ ); ///\endcond public: /// Calculate the temporal PSD for a Fourier mode in a multi-layer model. - /** + /** `PSD` and `freq` must have the same nonzero size. Frequencies must be finite, nonnegative, and strictly + * increasing. The AO system, integration controls, and complete atmosphere are validated before calculation. A + * precondition failure leaves `PSD` unchanged and clears `report` when supplied. * * \tparam parallel controls whether layers are calculated in parallel. Default is true. Set to false if this is * called inside a parallelized loop, as in \ref makePSDGrid. * - * \todo implement error checking - * \todo handle reports of convergence failures form singleLayerPSD when implemented. + * In permissive mode, recognized convergence failures are retained and summarized in `report`. Strict mode returns + * `error_t::liberr` if any layer has such a failure; the output PSD is incomplete and must be discarded whenever + * this function returns an error. * + * \returns `error_t::noerror` on success, or the first layer error in atmospheric-layer order. */ template - int multiLayerPSD( - std::vector &PSD, ///< [out] the calculated PSD - std::vector - &freq, ///< [in] the populated temporal frequency grid defining at which frequencies the PSD is calculated - realT m, ///< [in] the first index of the spatial frequency - realT n, ///< [in] the second index of the spatial frequency - int p, ///< [in] sets which mode is calculated (if basic modes, p = -1 for sine, p = +1 for cosine) - realT fmax = - 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The PSD is filled in - /// with a -17/3 power law past this frequency. If 0, then it is taken to be 150 Hz + 2*fastestPeak(m,n). - ); + error_t multiLayerPSD( + std::vector &PSD, /**< [out] calculated PSD */ + std::vector &freq, /**< [in] temporal-frequency grid */ + realT m, /**< [in] first spatial-frequency index */ + realT n, /**< [in] second spatial-frequency index */ + int p, /**< [in] Fourier-mode parity */ + realT fmax = 0, /**< [in] maximum exactly integrated frequency, or 0 for the default cutoff */ + reportT *report = nullptr, /**< [out] optional quadrature report */ + fourierTemporalPSDPolicy policy = fourierTemporalPSDPolicy::permissive /**< [in] non-convergence policy */ ); /// Calculate PSDs over a grid of spatial frequencies. /** The grid of spatial frequencies is square, set by the maximum value of m and n. @@ -317,14 +621,20 @@ struct fourierTemporalPSD * this adds overhead and cfitisio handles parallelization poorly due to the limitation on number of created file * pointers. * + * Inputs and AO-system state are validated before any output is created. A positive `fmax` switches each PSD to + * its asymptotic power-law tail above that frequency; zero selects the multilayer default cutoff. Calculation and + * write failures are collected by spatial-mode index and the first failure in grid order is returned after the + * parallel loop. Files completed before a calculation or write failure are retained. + * + * \returns `error_t::noerror` when the complete grid is written, or a typed validation, calculation, or output + * error. + * */ - void makePSDGrid( const std::string &dir, ///< [in] the directory for output of the PSDs - int mnMax, ///< [in] the maximum value of m and n in the grid. - realT dFreq, ///< [in] the temporal frequency spacing. - realT maxFreq, ///< [in] the maximum temporal frequency to calculate - realT fmax = 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The - ///< PSD is filled in with a -17/3 power law past - /// this frequency. If 0, then it is taken to be 150 Hz + 2*fastestPeak(m,n). + error_t makePSDGrid( const std::string &dir, ///< [in] the directory for output of the PSDs + int mnMax, ///< [in] the positive maximum value of m and n in the grid + realT dFreq, ///< [in] the positive temporal frequency spacing + realT maxFreq, ///< [in] the positive maximum temporal frequency to calculate + realT fmax = 0 ///< [in] maximum exactly calculated frequency, or 0 for the default cutoff ); /// Analyze a PSD grid under closed-loop control. @@ -413,36 +723,172 @@ fourierTemporalPSD::fourierTemporalPSD() } template -fourierTemporalPSD::fourierTemporalPSD( bool allocate ) +fourierTemporalPSD::fourierTemporalPSD( fourierTemporalPSD_detail::gslWorkspaceAllocator allocator ) + : m_workspaceAllocator( allocator ) { m_aosys = nullptr; initialize(); - - if( allocate ) - { - _w = gsl_integration_workspace_alloc( WSZ ); - } } template -fourierTemporalPSD::~fourierTemporalPSD() +error_t fourierTemporalPSD::allocateWorkspace() { - if( _w ) + if( m_workspace != nullptr ) + { + return error_t::noerror; + } + + if( m_workspaceAllocator == nullptr ) { - gsl_integration_workspace_free( _w ); + return internal::mxlib_error_report( error_t::invalidconfig, "GSL workspace allocator is null" ); } + + m_workspace.reset( m_workspaceAllocator( WSZ ) ); + if( m_workspace == nullptr ) + { + return internal::mxlib_error_report( error_t::allocerr, "could not allocate GSL integration workspace" ); + } + + return error_t::noerror; } template void fourierTemporalPSD::initialize() { _useBasis = basis::modified; - _w = 0; _absTol = 1e-10; _relTol = 1e-4; } +template +error_t fourierTemporalPSD::validateAtmosphere( int layer_i ) +{ + if( m_aosys == nullptr ) + { + return internal::mxlib_error_report( error_t::invalidconfig, "AO system pointer is null" ); + } + + if( !math::isFinite( m_aosys->D() ) || m_aosys->D() <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "AO system aperture diameter must be finite and positive" ); + } + + auto &atmosphere = m_aosys->atm; + const error_t atmosphereStatus = atmosphere.validate(); + if( atmosphereStatus != error_t::noerror ) + { + return atmosphereStatus; + } + + const size_t layerCount = atmosphere.n_layers(); + for( size_t index = 0; index < layerCount; ++index ) + { + if( atmosphere.layer_v_wind( static_cast( index ) ) <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "Fourier temporal PSD layers require positive wind speed" ); + } + } + + if( layer_i < -1 || ( layer_i >= 0 && static_cast( layer_i ) >= layerCount ) ) + { + return internal::mxlib_error_report( error_t::invalidarg, "atmosphere layer index is out of range" ); + } + + return error_t::noerror; +} + +template +error_t fourierTemporalPSD::validatePsdInputs( const std::vector &PSD, + const std::vector &freq, + realT m, + realT n, + int p, + realT fmax, + int layer_i, + fourierTemporalPSDPolicy policy ) +{ + if( freq.empty() || PSD.size() != freq.size() ) + { + return internal::mxlib_error_report( error_t::sizeerr, + "PSD and frequency vectors must have the same nonzero size" ); + } + + for( size_t index = 0; index < freq.size(); ++index ) + { + if( !math::isFinite( freq[index] ) || freq[index] < 0 || ( index > 0 && freq[index] <= freq[index - 1] ) ) + { + return internal::mxlib_error_report( + error_t::invalidarg, + "frequency grid must be finite, nonnegative, and strictly increasing" ); + } + } + + if( !math::isFinite( m ) || !math::isFinite( n ) || !math::isFinite( fmax ) || fmax < 0 ) + { + return internal::mxlib_error_report( + error_t::invalidarg, + "mode coordinates must be finite and frequency cutoff must be finite and nonnegative" ); + } + + if( p != -1 && p != 1 ) + { + return internal::mxlib_error_report( error_t::invalidarg, "Fourier-mode parity must be -1 or +1" ); + } + + if( _useBasis != basis::basic && _useBasis != basis::modified ) + { + return internal::mxlib_error_report( error_t::invalidarg, "value of _useBasis is not valid" ); + } + + if( policy != fourierTemporalPSDPolicy::permissive && policy != fourierTemporalPSDPolicy::strict ) + { + return internal::mxlib_error_report( error_t::invalidarg, "quadrature policy is not valid" ); + } + + if( !math::isFinite( _absTol ) || _absTol <= 0 || !math::isFinite( _relTol ) || _relTol <= 0 || _relTol >= 1 ) + { + return internal::mxlib_error_report( + error_t::invalidconfig, + "GSL absolute tolerance must be positive and relative tolerance must be between zero and one" ); + } + + if( !math::isFinite( m_f0 ) || m_f0 < 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "turbulence boiling parameter must be finite and nonnegative" ); + } + + const error_t atmosphereStatus = validateAtmosphere( layer_i ); + if( atmosphereStatus != error_t::noerror ) + { + return atmosphereStatus; + } + + if( _useBasis == basis::modified ) + { + if( !math::isFinite( m_aosys->lam_sci() ) || m_aosys->lam_sci() <= 0 || !math::isFinite( m_aosys->lam_wfs() ) || + m_aosys->lam_wfs() <= 0 || !math::isFinite( m_aosys->zeta() ) || + std::abs( m_aosys->zeta() ) >= math::half_pi() ) + { + return internal::mxlib_error_report( + error_t::invalidconfig, + "AO wavelengths must be positive and zenith angle must lie strictly between -pi/2 and pi/2" ); + } + } + + if( !math::isFinite( m_aosys->spatialFilter_ku() ) || m_aosys->spatialFilter_ku() <= 0 || + !math::isFinite( m_aosys->spatialFilter_kv() ) || m_aosys->spatialFilter_kv() <= 0 ) + { + return internal::mxlib_error_report( error_t::invalidconfig, + "AO spatial-filter limits must be finite and positive" ); + } + + return error_t::noerror; +} + template void fourierTemporalPSD::absTol( realT at ) { @@ -491,12 +937,51 @@ realT fourierTemporalPSD::fastestPeak( int m, int n ) } template -int fourierTemporalPSD::singleLayerPSD( - std::vector &PSD, std::vector &freq, realT m, realT n, int layer_i, int p, realT fmax ) +error_t fourierTemporalPSD::singleLayerPSD( std::vector &PSD, + std::vector &freq, + realT m, + realT n, + int layer_i, + int p, + realT fmax, + reportT *report, + fourierTemporalPSDPolicy policy ) +{ + reportT localReport; + reportT &activeReport = report == nullptr ? localReport : *report; + activeReport.clear(); + + const error_t status = validatePsdInputs( PSD, freq, m, n, p, fmax, layer_i, policy ); + if( status != error_t::noerror ) + { + return status; + } + + fourierTemporalPSD_detail::scopedGslErrorHandlerOff handlerGuard; + return singleLayerPSDImpl( PSD, freq, m, n, layer_i, p, fmax, activeReport, policy ); +} + +template +error_t fourierTemporalPSD::singleLayerPSDImpl( std::vector &PSD, + std::vector &freq, + realT m, + realT n, + int layer_i, + int p, + realT fmax, + reportT &report, + fourierTemporalPSDPolicy policy ) { if( fmax == 0 ) fmax = freq[freq.size() - 1]; + if( freq[0] > fmax ) + { + return internal::mxlib_error_report( + error_t::invalidarg, + "at least one exact frequency bin is required to initialize the PSD tail" ); + } + realT v_wind = m_aosys->atm.layer_v_wind( layer_i ); realT q_wind = m_aosys->atm.layer_dir( layer_i ); @@ -506,11 +991,13 @@ int fourierTemporalPSD::singleLayerPSD( realT scale = 2 * ( 1 / v_wind ); // Factor of 2 for negative frequencies. - // We'll get the occasional failure to reach tolerance error, just ignore them all for now. - gsl_set_error_handler_off(); - // Create a local instance so that we're reentrant - fourierTemporalPSD params( true ); + fourierTemporalPSD params( m_workspaceAllocator ); + const error_t workspaceStatus = params.allocateWorkspace(); + if( workspaceStatus != error_t::noerror ) + { + return workspaceStatus; + } params.m_aosys = m_aosys; params._layer_i = layer_i; @@ -529,21 +1016,22 @@ int fourierTemporalPSD::singleLayerPSD( params.m_modeCoeffs = m_modeCoeffs; params.m_minCoeffVal = m_minCoeffVal; - realT result, error; + realT result{ 0 }; + realT error{ 0 }; + error_t returnStatus = error_t::noerror; // Setup the GSL calculation gsl_function func; switch( _useBasis ) { - case basis::basic: // MXAO_FTPSD_BASIS_BASIC: - func.function = &F_basic; - break; - case basis::modified: // MXAO_FTPSD_BASIS_MODIFIED: - func.function = &F_mod; - break; - default: - internal::mxlib_error_report(error_t::invalidarg,"value of _useBasis is not valid." ); - return -1; + case basis::basic: // MXAO_FTPSD_BASIS_BASIC: + func.function = &F_basic; + break; + case basis::modified: // MXAO_FTPSD_BASIS_MODIFIED: + func.function = &F_mod; + break; + default: + return internal::mxlib_error_report( error_t::invalidarg, "value of _useBasis is not valid." ); } func.params = ¶ms; @@ -554,13 +1042,26 @@ int fourierTemporalPSD::singleLayerPSD( { params.m_f = freq[i]; - int ec = gsl_integration_qagi( &func, _absTol, _relTol, WSZ, params._w, &result, &error ); + const int ec = gsl_integration_qagi( &func, _absTol, _relTol, WSZ, params.m_workspace.get(), &result, &error ); + report.record( ec, static_cast( layer_i ), freq[i], result, error, _absTol, _relTol ); - if( ec == GSL_EDIVERGE ) + const error_t integrationStatus = fourierTemporalPSD_detail::applyPolicy( ec, policy ); + if( integrationStatus != error_t::noerror ) { - std::cerr << "GSL_EDIVERGE:" << p << " " << freq[i] << " " << v_wind << " " << m << " " << n << " " << m_m - << " " << m_n << "\n"; - std::cerr << "ignoring . . .\n"; + if( !fourierTemporalPSD_detail::isConvergenceStatus( ec ) ) + { + return internal::mxlib_error_report( integrationStatus, + std::string( "gsl_integration_qagi failed: " ) + + gsl_strerror( ec ) ); + } + + returnStatus = integrationStatus; + } + + if( !math::isFinite( result ) || !math::isFinite( error ) ) + { + return internal::mxlib_error_report( error_t::liberr, + "gsl_integration_qagi returned a nonfinite result or error estimate" ); } PSD[i] = scale * result; @@ -574,21 +1075,22 @@ int fourierTemporalPSD::singleLayerPSD( size_t j = i; if( j == freq.size() ) - return 0; + return returnStatus; - // First average result for last 50. - PSD[j] = - PSD[i - 50] * pow( freq[i - 50] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds()); - for( size_t k = 49; k > 0; --k ) + // First average up to the last 50 exactly integrated bins after projecting them to the first tail frequency. + constexpr size_t maximumTailAverageCount = 50; + const size_t tailAverageCount = std::min( i, maximumTailAverageCount ); + PSD[j] = 0; + for( size_t k = tailAverageCount; k > 0; --k ) { PSD[j] += PSD[i - k] * pow( freq[i - k] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds()); } - PSD[j] /= 50.0; + PSD[j] /= static_cast( tailAverageCount ); ++j; ++i; if( j == freq.size() ) - return 0; + return returnStatus; while( j < freq.size() ) { PSD[j] = @@ -596,15 +1098,26 @@ int fourierTemporalPSD::singleLayerPSD( ++j; } - return 0; + return returnStatus; } template -int fourierTemporalPSD::m_multiLayerPSD( - std::vector &PSD, std::vector &freq, realT m, realT n, int p, realT fmax, isParallel parallel ) +error_t fourierTemporalPSD::multiLayerPSDImpl( std::vector &PSD, + std::vector &freq, + realT m, + realT n, + int p, + realT fmax, + reportT &report, + fourierTemporalPSDPolicy policy, + isParallel parallel ) { static_cast( parallel ); + const size_t layerCount = m_aosys->atm.n_layers(); + std::vector layerStatus( layerCount, error_t::noerror ); + std::vector layerReport( layerCount ); + #pragma omp parallel { // Records each layer PSD @@ -613,48 +1126,98 @@ int fourierTemporalPSD::m_multiLayerPSD( #pragma omp for for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i ) { - singleLayerPSD( single_PSD, freq, m, n, i, p, fmax ); + std::fill( single_PSD.begin(), single_PSD.end(), 0 ); + layerStatus[i] = + singleLayerPSDImpl( single_PSD, freq, m, n, static_cast( i ), p, fmax, layerReport[i], policy ); // Now add the single layer PSD to the overall PSD, weighted by Cn2 #pragma omp critical - for( size_t j = 0; j < freq.size(); ++j ) + if( layerStatus[i] == error_t::noerror ) { - PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j]; + for( size_t j = 0; j < freq.size(); ++j ) + { + PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j]; + } } } } - return 0; + error_t returnStatus = error_t::noerror; + for( size_t i = 0; i < layerCount; ++i ) + { + report.merge( layerReport[i] ); + if( returnStatus == error_t::noerror && layerStatus[i] != error_t::noerror ) + { + returnStatus = layerStatus[i]; + } + } + + return returnStatus; } template -int fourierTemporalPSD::m_multiLayerPSD( - std::vector &PSD, std::vector &freq, realT m, realT n, int p, realT fmax, isParallel parallel ) +error_t fourierTemporalPSD::multiLayerPSDImpl( std::vector &PSD, + std::vector &freq, + realT m, + realT n, + int p, + realT fmax, + reportT &report, + fourierTemporalPSDPolicy policy, + isParallel parallel ) { static_cast( parallel ); // Records each layer PSD std::vector single_PSD( freq.size() ); + error_t returnStatus = error_t::noerror; for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i ) { - singleLayerPSD( single_PSD, freq, m, n, i, p, fmax ); + std::fill( single_PSD.begin(), single_PSD.end(), 0 ); + reportT layerReport; + const error_t layerStatus = + singleLayerPSDImpl( single_PSD, freq, m, n, static_cast( i ), p, fmax, layerReport, policy ); + report.merge( layerReport ); + if( returnStatus == error_t::noerror && layerStatus != error_t::noerror ) + { + returnStatus = layerStatus; + } // Now add the single layer PSD to the overall PSD, weighted by Cn2 - for( size_t j = 0; j < freq.size(); ++j ) + if( layerStatus == error_t::noerror ) { - PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j]; + for( size_t j = 0; j < freq.size(); ++j ) + { + PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j]; + } } } - return 0; + return returnStatus; } template template -int fourierTemporalPSD::multiLayerPSD( - std::vector &PSD, std::vector &freq, realT m, realT n, int p, realT fmax ) +error_t fourierTemporalPSD::multiLayerPSD( std::vector &PSD, + std::vector &freq, + realT m, + realT n, + int p, + realT fmax, + reportT *report, + fourierTemporalPSDPolicy policy ) { + reportT localReport; + reportT &activeReport = report == nullptr ? localReport : *report; + activeReport.clear(); + + const error_t validationStatus = validatePsdInputs( PSD, freq, m, n, p, fmax, -1, policy ); + if( validationStatus != error_t::noerror ) + { + return validationStatus; + } + // PSD is zeroed every time to make sure we don't accumulate on repeated calls for( size_t j = 0; j < PSD.size(); ++j ) PSD[j] = 0; @@ -664,13 +1227,48 @@ int fourierTemporalPSD::multiLayerPSD( fmax = 150 + 2 * fastestPeak( m, n ); } - return m_multiLayerPSD( PSD, freq, m, n, p, fmax, isParallel() ); + fourierTemporalPSD_detail::scopedGslErrorHandlerOff handlerGuard; + return multiLayerPSDImpl( PSD, freq, m, n, p, fmax, activeReport, policy, isParallel() ); } template -void fourierTemporalPSD::makePSDGrid( +error_t fourierTemporalPSD::makePSDGrid( const std::string &dir, int mnMax, realT dFreq, realT maxFreq, realT fmax ) { + if( dir.empty() ) + { + return internal::mxlib_error_report( error_t::invalidarg, "PSD grid output directory must not be empty" ); + } + + if( mnMax <= 0 || !math::isFinite( dFreq ) || dFreq <= 0 || !math::isFinite( maxFreq ) || maxFreq <= 0 || + !math::isFinite( fmax ) || fmax < 0 ) + { + return internal::mxlib_error_report( + error_t::invalidarg, + "PSD grid extent and frequency controls must be finite and positive, with a nonnegative cutoff" ); + } + + const realT sampleCount = maxFreq / dFreq; + if( !math::isFinite( sampleCount ) || sampleCount > static_cast( std::numeric_limits::max() ) ) + { + return internal::mxlib_error_report( error_t::sizeerr, "PSD grid sample count exceeds the supported range" ); + } + + const std::vector validationFrequency{ 0 }; + const std::vector validationPsd{ 0 }; + const error_t validationStatus = validatePsdInputs( validationPsd, + validationFrequency, + 0, + 0, + 1, + fmax, + -1, + fourierTemporalPSDPolicy::permissive ); + if( validationStatus != error_t::noerror ) + { + return validationStatus; + } + std::vector freq; std::vector spf; @@ -685,11 +1283,19 @@ void fourierTemporalPSD::makePSDGrid( N += 1; /*** Dump Params to file ***/ - mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); + error_t status = ioutils::createDirectories( dir ); + if( status != error_t::noerror ) + { + return internal::mxlib_error_report( status, "could not create PSD grid output directory" ); + } std::ofstream fout; fn = dir + '/' + "params.txt"; fout.open( fn ); + if( !fout.is_open() ) + { + return internal::mxlib_error_report( error_t::fileoerr, "could not open PSD grid parameter file" ); + } fout << "#---------------------------\n"; m_aosys->dumpAOSystem( fout ); @@ -706,19 +1312,31 @@ void fourierTemporalPSD::makePSDGrid( fout << "#---------------------------\n"; fout.close(); + if( !fout ) + { + return internal::mxlib_error_report( error_t::filewerr, "could not write PSD grid parameter file" ); + } // Make directory std::string psddir = dir + "/psds"; - mkdir( psddir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); + status = ioutils::createDirectories( psddir ); + if( status != error_t::noerror ) + { + return internal::mxlib_error_report( status, "could not create PSD output directory" ); + } // Create frequency scale. math::vectorScale( freq, N, dFreq, 0 ); // dFreq); //offset from 0 by dFreq, so f=0 not included fn = psddir + '/' + "freq.binv"; - ioutils::writeBinVector( fn, freq ); + if( ioutils::writeBinVector( fn, freq ) != 0 ) + { + return internal::mxlib_error_report( error_t::filewerr, "could not write PSD frequency grid" ); + } size_t nLoops = 0.5 * spf.size(); + std::vector modeStatus( nLoops, error_t::noerror ); ipc::ompLoopWatcher<> watcher( nLoops, std::cout ); @@ -743,16 +1361,35 @@ void fourierTemporalPSD::makePSDGrid( continue; } - multiLayerPSD( PSD, freq, m, n, 1, fmax ); + modeStatus[i] = multiLayerPSD( PSD, freq, m, n, 1, fmax ); + if( modeStatus[i] != error_t::noerror ) + { + watcher.incrementAndOutputStatus(); + continue; + } - fname = std::format("{}/psd_{}_{}.binv",psddir,m,n); - // psddir + '/' + "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + fname = std::format( "{}/psd_{}_{}.binv", psddir, m, n ); + // psddir + '/' + "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + + // ".binv"; - ioutils::writeBinVector( fname, PSD ); + if( ioutils::writeBinVector( fname, PSD ) != 0 ) + { + modeStatus[i] = error_t::filewerr; + } watcher.incrementAndOutputStatus(); } } + + for( size_t index = 0; index < modeStatus.size(); ++index ) + { + if( modeStatus[index] != error_t::noerror ) + { + return modeStatus[index]; + } + } + + return error_t::noerror; } template @@ -843,14 +1480,14 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir { for( size_t s = 0; s < mags.size(); ++s ) { - std::string psdOutDir = std::format("{}/outputPSDS_{}_si",dir, mags[s]); - //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si"; + std::string psdOutDir = std::format( "{}/outputPSDS_{}_si", dir, mags[s] ); + // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si"; mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); if( doLP ) { - std::string psdOutDir = std::format("{}/outputPSDS_{}_lp",dir, mags[s]); - //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp"; + std::string psdOutDir = std::format( "{}/outputPSDS_{}_lp", dir, mags[s] ); + // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp"; mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); } } @@ -859,6 +1496,7 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir m_aosys->beta_p( 1, 1 ); ipc::ompLoopWatcher<> watcher( nModes * mags.size(), std::cout ); + std::atomic analysisStatus{ static_cast( error_t::noerror ) }; for( size_t s = 0; s < mags.size(); ++s ) { @@ -976,8 +1614,8 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir if( writeXfer ) { - std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; + std::string tfOutFile = std::format( "{}/outputTF_{}_si", dir, mags[s] ); + // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; ioutils::createDirectories( tfOutFile ); } @@ -985,8 +1623,8 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir { if( writeXfer ) { - std::string tfOutFile = std::format("{}/outputTF_{}_lp", dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; + std::string tfOutFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] ); + // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; ioutils::createDirectories( tfOutFile ); } } @@ -1103,7 +1741,14 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; else { // Calculate gain using the POL PSD - gopt = go_si.optGainOpenLoop( var, tPSDpPOL, tPSDn, gmax, true ); + error_t gainStatus = go_si.optGainOpenLoop( gopt, var, tPSDpPOL, tPSDn, true ); + if( gainStatus != error_t::noerror ) + { + int expected = static_cast( error_t::noerror ); + analysisStatus.compare_exchange_strong( expected, static_cast( gainStatus ) ); + gopt = 0; + var = go_si.clVariance( tPSDp, tPSDn, gopt ); + } if( m_uncorrectedOG ) { @@ -1139,21 +1784,19 @@ 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 << __FILE__ << ' ' << __LINE__ << '\n'; + int expected = static_cast( error_t::noerror ); + analysisStatus.compare_exchange_strong( expected, static_cast( rv ) ); } for( int n = 0; n < lpNc; ++n ) @@ -1199,7 +1842,7 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( gopt_lp > gopt && var_lp > var ) { - //Set LP to SI (or off if SI is off) + // Set LP to SI (or off if SI is off) gopt_lp = gopt; var_lp = var; go_lp.a( std::vector( { 1 } ) ); @@ -1245,17 +1888,17 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( writeXfer ) { - std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]); + std::string tfOutFile = std::format( "{}/outputTF_{}_si", dir, mags[s] ); // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; - std::string etfOutFile = std::format("{}/etf_{}_{}.binv", tfOutFile,m, n); - //tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' + - // ioutils::convert ToString( n ) + ".binv"; + std::string etfOutFile = std::format( "{}/etf_{}_{}.binv", tfOutFile, m, n ); + // tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' + + // ioutils::convert ToString( n ) + ".binv"; ioutils::writeBinVector( etfOutFile, ETFxn ); - std::string ntfOutFile = std::format("{}/ntf_{}_{}.binv",tfOutFile, m, n); - //tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' + - // ioutils::convert ToString( n ) + ".binv"; + std::string ntfOutFile = std::format( "{}/ntf_{}_{}.binv", tfOutFile, m, n ); + // tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' + + // ioutils::convert ToString( n ) + ".binv"; ioutils::writeBinVector( ntfOutFile, NTFxn ); if( i == 0 ) // Write freq on the first one @@ -1300,17 +1943,17 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( writeXfer ) { - std::string tfOutFile = std::format("{}/outputTF_{}_lp",dir,mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; + std::string tfOutFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] ); + // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; - std::string etfOutFile = std::format("{}/etf_{}_{}.binv", tfOutFile, m, n); - //tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' + - // ioutils::convert ToString( n ) + ".binv"; + std::string etfOutFile = std::format( "{}/etf_{}_{}.binv", tfOutFile, m, n ); + // tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' + + // ioutils::convert ToString( n ) + ".binv"; ioutils::writeBinVector( etfOutFile, ETFxn ); - std::string ntfOutFile = std::format("{}/ntf_{}_{}.binv", tfOutFile, m, n); - //tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' + - // ioutils::convert ToString( n ) + ".binv"; + std::string ntfOutFile = std::format( "{}/ntf_{}_{}.binv", tfOutFile, m, n ); + // tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' + + // ioutils::convert ToString( n ) + ".binv"; ioutils::writeBinVector( ntfOutFile, NTFxn ); if( i == 0 ) // Write freq on the first one @@ -1340,11 +1983,12 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; // Calculate the controlled PSDs and output if( writePSDs ) { - std::string psdOutFile = std::format("{}/outputPSDs_{}_si/psd_{}_{}.binv",dir,mags[s],m,n); + std::string psdOutFile = + std::format( "{}/outputPSDs_{}_si/psd_{}_{}.binv", dir, mags[s], m, n ); - // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/"; - //psdOutFile += - // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/"; + // psdOutFile += + // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; std::vector psdOut( tPSDp.size() + tPSDpHF.size() ); @@ -1381,18 +2025,20 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( i == 0 ) // Write freq on the first one { - psdOutFile = std::format("{}/outputPSDs_{}_si/freq.binv", dir, mags[s]); + psdOutFile = std::format( "{}/outputPSDs_{}_si/freq.binv", dir, mags[s] ); // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/freq.binv"; - //ioutils::writeBinVector( psdOutFile, tfreqHF ); + // ioutils::writeBinVector( psdOutFile, tfreqHF ); } if( doLP ) { - std::string psdOutFile = std::format("{}/outputPSDs_{}_lp/psd_{}_{}.binv",dir,mags[s],m,n); + std::string psdOutFile = + std::format( "{}/outputPSDs_{}_lp/psd_{}_{}.binv", dir, mags[s], m, n ); - //psdOutFile = dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/"; - //psdOutFile += - // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + // psdOutFile = dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/"; + // psdOutFile += + // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + + // ".binv"; // Calculate the output PSD if gains are applied if( gopt_lp > 0 ) @@ -1426,10 +2072,11 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( i == 0 ) { - psdOutFile = std::format("{}/outputPSDs_{}_lp/freq.binv", dir, mags[s]); - //psdOutFile = - // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/freq.binv"; - //ioutils::writeBinVector( psdOutFile, tfreq ); + psdOutFile = std::format( "{}/outputPSDs_{}_lp/freq.binv", dir, mags[s] ); + // psdOutFile = + // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + + // "_lp/freq.binv"; + // ioutils::writeBinVector( psdOutFile, tfreq ); } } } @@ -1439,15 +2086,20 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; } // omp for i..nModes } // omp Parallel + if( analysisStatus.load() != static_cast( error_t::noerror ) ) + { + return analysisStatus.load(); + } + Eigen::Array cim; fits::fitsFile ff; - std::string fn = std::format("{}/gainmap_{}_si.fits",dir, mags[s]); - //dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + std::string fn = std::format( "{}/gainmap_{}_si.fits", dir, mags[s] ); + // dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, gains ); - fn = std::format("{}/varmap_{}_si.fits",dir, mags[s]); - //dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + fn = std::format( "{}/varmap_{}_si.fits", dir, mags[s] ); + // dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, vars ); cim = vars; @@ -1456,29 +2108,29 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; S_si.push_back( strehl ); cim /= strehl; - fn = std::format("{}/contrast_{}_si.fits",dir, mags[s]); - //dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + fn = std::format( "{}/contrast_{}_si.fits", dir, mags[s] ); + // dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, cim ); if( lifetimeTrials > 0 ) { - fn = std::format("{}/speckleLifetimes_{}_si.fits",dir, mags[s]); - //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + fn = std::format( "{}/speckleLifetimes_{}_si.fits", dir, mags[s] ); + // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, speckleLifetimes ); } if( doLP ) { - fn = std::format("{}/gainmap_{}_lp.fits",dir, mags[s]); - //dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/gainmap_{}_lp.fits", dir, mags[s] ); + // dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, gains_lp ); - fn = std::format("{}/lpcmap_{}_lp.fits",dir, mags[s]); - //dir + "/lpcmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/lpcmap_{}_lp.fits", dir, mags[s] ); + // dir + "/lpcmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, lpC ); - fn = std::format("{}/varmap_{}_lp.fits",dir, mags[s]); - //dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/varmap_{}_lp.fits", dir, mags[s] ); + // dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, vars_lp ); cim = vars_lp; @@ -1489,14 +2141,14 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; S_lp.push_back( Slp ); cim /= Slp; - fn = std::format("{}/contrast_{}_lp.fits",dir, mags[s]); - //dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/contrast_{}_lp.fits", dir, mags[s] ); + // dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, cim ); if( lifetimeTrials > 0 ) { - fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]); - //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] ); + // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, speckleLifetimes_lp ); } } @@ -1744,35 +2396,35 @@ int fourierTemporalPSD::intensityPSD( if( inside ) { - tfInFile = std::format("{}/outputTF_{}_si",dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; + tfInFile = std::format( "{}/outputTF_{}_si", dir, mags[s] ); + // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; - etfInFile = std::format("{}/etf_{}_{}.binv",tfInFile, m, n); - //tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + etfInFile = std::format( "{}/etf_{}_{}.binv", tfInFile, m, n ); + // tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; ioutils::readBinVector( tPSDc, etfInFile ); sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order for( size_t j = 0; j < psd2sidedc.size(); ++j ) ETFsi[i][j] = psd2sidedc[j]; - ntfInFile = std::format("{}/ntf_{}_{}.binv",tfInFile, m, n); - //tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + ntfInFile = std::format( "{}/ntf_{}_{}.binv", tfInFile, m, n ); + // tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; ioutils::readBinVector( tPSDc, ntfInFile ); sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order for( size_t j = 0; j < psd2sidedc.size(); ++j ) NTFsi[i][j] = psd2sidedc[j]; - tfInFile = std::format("{}/outputTF_{}_lp",dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; + tfInFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] ); + // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; - etfInFile = std::format("{}/etf_{}_{}.binv",tfInFile, m, n); - //tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + etfInFile = std::format( "{}/etf_{}_{}.binv", tfInFile, m, n ); + // tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; ioutils::readBinVector( tPSDc, etfInFile ); sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order for( size_t j = 0; j < psd2sidedc.size(); ++j ) ETFlp[i][j] = psd2sidedc[j]; - ntfInFile = std::format("{}/ntf_{}_{}.binv",tfInFile, m, n); - //tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + ntfInFile = std::format( "{}/ntf_{}_{}.binv", tfInFile, m, n ); + // tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; ioutils::readBinVector( tPSDc, ntfInFile ); sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order for( size_t j = 0; j < psd2sidedc.size(); ++j ) @@ -2195,22 +2847,22 @@ int fourierTemporalPSD::intensityPSD( /*********************************************************************/ // 4.0) Write the results to disk /*********************************************************************/ - fn = std::format("{}/speckleLifetimes_{}_si.fits",dir, mags[s]); - //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + fn = std::format( "{}/speckleLifetimes_{}_si.fits", dir, mags[s] ); + // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, taus ); - fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]); - //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] ); + // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, tauslp ); if( writePSDs ) { - fn = std::format("{}/specklePSDs_{}_si.fits",dir, mags[s]); - //dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; + fn = std::format( "{}/specklePSDs_{}_si.fits", dir, mags[s] ); + // dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_si.fits"; ff.write( fn, imc ); - fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]); - //dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; + fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] ); + // dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_lp.fits"; ff.write( fn, imclp ); } @@ -2231,8 +2883,8 @@ template int fourierTemporalPSD::getGridPSD( std::vector &psd, const std::string &dir, int m, int n ) { std::string fn; - fn = std::format("{}/psds/psd_{}_{}.binv",dir,m,n); - //dir + "/psds/psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; + fn = std::format( "{}/psds/psd_{}_{}.binv", dir, m, n ); + // dir + "/psds/psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv"; return ioutils::readBinVector( psd, fn ); } @@ -2421,11 +3073,10 @@ realT F_mod( realT kv, void *params ) P2 *= QQ; - return 0.5*(P1 + P2); + return 0.5 * ( P1 + P2 ); } } - /*extern template struct fourierTemporalPSD, std::ostream>>;*/ 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/source/ao/analysis/clGainOpt.cpp b/source/ao/analysis/clGainOpt.cpp index ffffa829f..593b92ea4 100644 --- a/source/ao/analysis/clGainOpt.cpp +++ b/source/ao/analysis/clGainOpt.cpp @@ -40,90 +40,108 @@ namespace impl { template -realT _optGainOpenLoop( clGainOptOptGain_OL &olgo, - realT &var, - const realT &gmax, - const realT &minFindMin, - const realT &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ) +mx::error_t _optGainOpenLoop( realT &gain, + realT &var, + clGainOptOptGain_OL &olgo, + const realT &minimumGain, + const realT &maximumGain, + int minFindBits, + uintmax_t minFindMaxIter, + uintmax_t &iters ) { - realT gopt; - + gain = std::numeric_limits::quiet_NaN(); + var = std::numeric_limits::quiet_NaN(); iters = minFindMaxIter; try { std::pair brack; - brack = boost::math::tools::brent_find_minima, realT>( - olgo, minFindMin, minFindMaxFact * gmax, minFindBits, iters ); - gopt = brack.first; + brack = boost::math::tools::brent_find_minima, realT>( olgo, + minimumGain, + maximumGain, + minFindBits, + iters ); + gain = brack.first; var = brack.second; } catch( ... ) { - std::cerr << "optGainOpenLoop: No root found\n"; - gopt = minFindMaxFact * gmax; - var = 0; + return error_t::exception; + } + + if( iters >= minFindMaxIter ) + { + return error_t::timeout; } - return gopt; + return error_t::noerror; } template <> -float optGainOpenLoop( clGainOptOptGain_OL &olgo, - float &var, - const float &gmax, - const float &minFindMin, - const float &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ) +mx::error_t optGainOpenLoop( float &gain, + float &var, + clGainOptOptGain_OL &olgo, + const float &minimumGain, + const float &maximumGain, + int minFindBits, + uintmax_t minFindMaxIter, + uintmax_t &iters ) { - return _optGainOpenLoop( olgo, var, gmax, minFindMin, minFindMaxFact, minFindBits, minFindMaxIter, iters ); + return _optGainOpenLoop( gain, var, olgo, minimumGain, maximumGain, minFindBits, minFindMaxIter, iters ); } template <> -double optGainOpenLoop( clGainOptOptGain_OL &olgo, - double &var, - const double &gmax, - const double &minFindMin, - const double &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t &iters ) +mx::error_t optGainOpenLoop( double &gain, + double &var, + clGainOptOptGain_OL &olgo, + const double &minimumGain, + const double &maximumGain, + int minFindBits, + uintmax_t minFindMaxIter, + uintmax_t &iters ) { - return _optGainOpenLoop( olgo, var, gmax, minFindMin, minFindMaxFact, minFindBits, minFindMaxIter, iters ); + return _optGainOpenLoop( gain, var, olgo, minimumGain, maximumGain, minFindBits, minFindMaxIter, iters ); } template <> -long double optGainOpenLoop( clGainOptOptGain_OL &olgo, +mx::error_t optGainOpenLoop( long double &gain, long double &var, - const long double &gmax, - const long double &minFindMin, - const long double &minFindMaxFact, + clGainOptOptGain_OL &olgo, + const long double &minimumGain, + const long double &maximumGain, int minFindBits, uintmax_t minFindMaxIter, uintmax_t &iters ) { - return _optGainOpenLoop( - olgo, var, gmax, minFindMin, minFindMaxFact, minFindBits, minFindMaxIter, iters ); + return _optGainOpenLoop( gain, + var, + olgo, + minimumGain, + maximumGain, + minFindBits, + minFindMaxIter, + iters ); } #ifdef HASQUAD template <> -__float128 optGainOpenLoop<__float128>( clGainOptOptGain_OL<__float128> &olgo, - __float128 &var, - __float128 &gmax, - __float128 &minFindMin, - __float128 &minFindMaxFact, - int minFindBits, - uintmax_t minFindMaxIter, - uintmax_t iters ) +mx::error_t optGainOpenLoop<__float128>( __float128 &gain, + __float128 &var, + clGainOptOptGain_OL<__float128> &olgo, + const __float128 &minimumGain, + const __float128 &maximumGain, + int minFindBits, + uintmax_t minFindMaxIter, + uintmax_t &iters ) { - return _optGainOpenLoop<__float128>( - olgo, var, gmax, minFindMin, minFindMaxFact, minFindBits, minFindMaxIter, iters ); + return _optGainOpenLoop<__float128>( gain, + var, + olgo, + minimumGain, + maximumGain, + minFindBits, + minFindMaxIter, + iters ); } #endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fc6ba464b..7444e6b3e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,9 @@ 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/clGainOpt_test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ao/analysis/fourierTemporalPSD_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 +16,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..005db08fa 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -6,6 +6,9 @@ INCLUDES += -I../include OBJS = testMain.o \ include/ao/analysis/aoAtmosphere_test.o \ + include/ao/analysis/clAOLinearPredictor_test.o \ + include/ao/analysis/clGainOpt_test.o \ + include/ao/analysis/fourierTemporalPSD_test.o \ include/ao/analysis/aoSystem_test.o \ include/ao/analysis/aoPSDs_test.o \ include/astro/astroDynamics_test.o \ @@ -13,6 +16,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/aoAtmosphere_test.cpp b/tests/include/ao/analysis/aoAtmosphere_test.cpp index 8c6800aa8..e9bd9aa6e 100644 --- a/tests/include/ao/analysis/aoAtmosphere_test.cpp +++ b/tests/include/ao/analysis/aoAtmosphere_test.cpp @@ -6,22 +6,22 @@ #include "../../../../include/ao/analysis/aoAtmosphere.hpp" +#include + typedef double realT; using namespace mx::app; using namespace mx::AO::analysis; -/** \test Scenario: Loading aoAtmosphere config settings - * - * Verify parsing of config - * - * \anchor tests_ao_analysis_aoAtmosphere_config - */ +/// Verify parsing and validation of atmosphere configuration settings. +/** Exercises mx::AO::analysis::aoAtmosphere::setupConfig and + * mx::AO::analysis::aoAtmosphere::loadConfig. */ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" ) { - GIVEN( "a valid config vile" ) + GIVEN( "a valid config file" ) { - aoAtmosphere atm; // This will be cumulative + aoAtmosphere atm; + atm.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); WHEN( "all normal settings" ) { @@ -44,7 +44,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.r_0() == 0.25 ); REQUIRE( atm.lam_0() == 0.4e-9 ); @@ -84,6 +84,8 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" } WHEN( "setting v_wind and z_mean" ) { + atm.L_0( { 25.0, 25.0, 25.0 } ); + atm.l_0( { 0.0, 0.0, 0.0 } ); appConfigurator config; writeConfigFile( "aoAtmosphere.conf", @@ -93,7 +95,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.v_wind() == Approx( 15.0 ) ); REQUIRE( atm.z_mean() == Approx( 5001.0 ) ); @@ -108,11 +110,11 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.nonKolmogorov() == false ); } - WHEN( "setting nonKolmogorov to true" ) + WHEN( "setting nonKolmogorov to false" ) { appConfigurator config; @@ -120,7 +122,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.nonKolmogorov() == true ); } @@ -136,7 +138,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.alpha( 0 ) == Approx( 4.7 ) ); REQUIRE( atm.nonKolmogorov() == true ); @@ -153,7 +155,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.beta( 0 ) == Approx( 0.026 ) ); REQUIRE( atm.nonKolmogorov() == true ); @@ -170,7 +172,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" atm.setupConfig( config ); config.readConfig( "aoAtmosphere.conf" ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.beta_0( 0 ) == Approx( 1e-7 ) ); REQUIRE( atm.nonKolmogorov() == true ); @@ -178,7 +180,8 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" } GIVEN( "command line options" ) { - aoAtmosphere atm; // This will be cumulative + aoAtmosphere atm; + atm.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); WHEN( "setting nonKolmogorov to true" ) { @@ -195,7 +198,7 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" REQUIRE( atm.nonKolmogorov() == false ); atm.setupConfig( config ); config.parseCommandLine( argvs.size(), argv ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.nonKolmogorov() == true ); } @@ -215,8 +218,207 @@ SCENARIO( "Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]" REQUIRE( atm.nonKolmogorov() == true ); atm.setupConfig( config ); config.parseCommandLine( argvs.size(), argv ); - atm.loadConfig( config ); + REQUIRE( atm.loadConfig( config ) == mx::error_t::noerror ); REQUIRE( atm.nonKolmogorov() == false ); } } } + +/// Verify that layer-strength normalization rejects invalid inputs without changing state. +/** Exercises mx::AO::analysis::aoAtmosphere::layer_Cn2. */ +TEST_CASE( "Atmosphere layer-strength normalization is transactional", "[ao::analysis::aoAtmosphere]" ) +{ + aoAtmosphere atmosphere; + atmosphere.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + const std::vector originalStrength = atmosphere.layer_Cn2(); + + SECTION( "empty strengths" ) + { + REQUIRE( atmosphere.layer_Cn2( std::vector{} ) == mx::error_t::invalidarg ); + } + + SECTION( "negative strength" ) + { + REQUIRE( atmosphere.layer_Cn2( { 1.0, -1.0 } ) == mx::error_t::invalidarg ); + } + + SECTION( "nonfinite strength" ) + { + REQUIRE( atmosphere.layer_Cn2( { 1.0, std::numeric_limits::quiet_NaN() } ) == mx::error_t::invalidarg ); + } + + SECTION( "nonfinite sum" ) + { + REQUIRE( atmosphere.layer_Cn2( { std::numeric_limits::max(), std::numeric_limits::max() } ) == + mx::error_t::invalidarg ); + } + + SECTION( "zero sum" ) + { + REQUIRE( atmosphere.layer_Cn2( { 0.0, 0.0 } ) == mx::error_t::invalidarg ); + } + + SECTION( "invalid reference wavelength" ) + { + REQUIRE( atmosphere.layer_Cn2( { 1.0 }, -1.0 ) == mx::error_t::invalidarg ); + } + + REQUIRE( atmosphere.layer_Cn2() == originalStrength ); + REQUIRE( atmosphere.validate() == mx::error_t::noerror ); +} + +/// Verify that complete atmosphere validation rejects every layer-vector mismatch. +/** Exercises mx::AO::analysis::aoAtmosphere::validate with each core and non-Kolmogorov vector setter. */ +TEST_CASE( "Atmosphere validation rejects mismatched layer vectors", "[ao::analysis::aoAtmosphere]" ) +{ + aoAtmosphere atmosphere; + atmosphere.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + SECTION( "layer strengths" ) + { + REQUIRE( atmosphere.layer_Cn2( { 1.0, 1.0 } ) == mx::error_t::noerror ); + } + + SECTION( "outer scales" ) + { + atmosphere.L_0( { 25.0, 25.0 } ); + } + + SECTION( "inner scales" ) + { + atmosphere.l_0( { 0.0, 0.0 } ); + } + + SECTION( "layer heights" ) + { + atmosphere.layer_z( { 0.0, 0.0 } ); + } + + SECTION( "wind speeds" ) + { + atmosphere.layer_v_wind( { 10.0, 10.0 } ); + } + + SECTION( "wind directions" ) + { + atmosphere.layer_dir( { 0.0, 0.0 } ); + } + + SECTION( "non-Kolmogorov exponents" ) + { + atmosphere.alpha( { 11.0 / 3.0, 11.0 / 3.0 } ); + } + + SECTION( "non-Kolmogorov normalizations" ) + { + atmosphere.beta( { 1.0, 1.0 } ); + } + + SECTION( "non-Kolmogorov constants" ) + { + atmosphere.beta_0( { 0.0, 0.0 } ); + } + + REQUIRE( atmosphere.validate() == mx::error_t::sizeerr ); +} + +/// Verify that atmosphere validation rejects nonphysical scalar and layer values. +/** Exercises mx::AO::analysis::aoAtmosphere::validate after mutation through its public setters. */ +TEST_CASE( "Atmosphere validation rejects invalid physical values", "[ao::analysis::aoAtmosphere]" ) +{ + aoAtmosphere atmosphere; + atmosphere.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + SECTION( "invalid Fried parameter" ) + { + atmosphere.r_0( 0.0, 0.5e-6 ); + } + + SECTION( "invalid observatory height" ) + { + atmosphere.h_obs( -1.0 ); + } + + SECTION( "invalid atmospheric scale height" ) + { + atmosphere.H( 0.0 ); + } + + SECTION( "nonfinite outer scale" ) + { + atmosphere.L_0( std::vector{ std::numeric_limits::infinity() } ); + } + + SECTION( "negative inner scale" ) + { + atmosphere.l_0( std::vector{ -1.0 } ); + } + + SECTION( "negative layer height" ) + { + atmosphere.layer_z( std::vector{ -1.0 } ); + } + + SECTION( "negative wind speed" ) + { + atmosphere.layer_v_wind( std::vector{ -1.0 } ); + } + + SECTION( "nonfinite wind direction" ) + { + atmosphere.layer_dir( std::vector{ std::numeric_limits::quiet_NaN() } ); + } + + SECTION( "invalid non-Kolmogorov normalization" ) + { + atmosphere.nonKolmogorov( true ); + atmosphere.beta( std::vector{ 0.0 } ); + } + + REQUIRE( atmosphere.validate() == mx::error_t::invalidconfig ); +} + +/// Verify that atmosphere validation preserves model-specific physical conventions. +/** Exercises mx::AO::analysis::aoAtmosphere::validate for static layers, infinite outer scale, and explicit + * non-Kolmogorov normalization. */ +TEST_CASE( "Atmosphere validation accepts supported physical conventions", "[ao::analysis::aoAtmosphere]" ) +{ + aoAtmosphere atmosphere; + atmosphere.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + SECTION( "static layer" ) + { + atmosphere.layer_v_wind( std::vector{ 0.0 } ); + } + + SECTION( "infinite outer scale" ) + { + atmosphere.L_0( std::vector{ 0.0 } ); + } + + SECTION( "non-Kolmogorov model without a Fried parameter" ) + { + atmosphere.r_0( 0.0, 0.5e-6 ); + atmosphere.nonKolmogorov( true ); + } + + REQUIRE( atmosphere.validate() == mx::error_t::noerror ); +} + +/// Verify that configuration loading returns layer-strength validation failures. +/** Exercises mx::AO::analysis::aoAtmosphere::loadConfig and + * mx::AO::analysis::aoAtmosphere::layer_Cn2. */ +TEST_CASE( "Atmosphere configuration propagates invalid layer strengths", "[ao::analysis::aoAtmosphere]" ) +{ + aoAtmosphere atmosphere; + atmosphere.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + appConfigurator config; + writeConfigFile( "aoAtmosphere.conf", { "atm" }, { "layer_Cn2" }, { "0,0" } ); + atmosphere.setupConfig( config ); + config.readConfig( "aoAtmosphere.conf" ); + + REQUIRE( atmosphere.loadConfig( config ) == mx::error_t::invalidarg ); + REQUIRE( atmosphere.validate() == mx::error_t::noerror ); + REQUIRE( atmosphere.layer_Cn2() == std::vector{ 1.0 } ); +} diff --git a/tests/include/ao/analysis/aoSystem_test.cpp b/tests/include/ao/analysis/aoSystem_test.cpp index a4afe3725..8a0793c00 100644 --- a/tests/include/ao/analysis/aoSystem_test.cpp +++ b/tests/include/ao/analysis/aoSystem_test.cpp @@ -11,11 +11,8 @@ typedef double realT; using namespace mx::app; using namespace mx::AO::analysis; -/** Scenario: Loading aoSystem config settings - * - * Verify parsing of config - * \anchor tests_ao_analysis_aoSystem_config - */ +/// Verify parsing of AO-system configuration settings. +/** Exercises mx::AO::analysis::aoSystem::setupConfig and mx::AO::analysis::aoSystem::loadConfig. */ SCENARIO( "Loading aoSystem config settings", "[ao::analysis::aoSystem]" ) { GIVEN( "no config file" ) @@ -210,3 +207,23 @@ SCENARIO( "Loading aoSystem config settings", "[ao::analysis::aoSystem]" ) } } } + +/// Verify that AO-system configuration loading propagates atmosphere validity. +/** Exercises mx::AO::analysis::aoSystem::loadConfig and mx::AO::analysis::aoAtmosphere::setSingleLayer. */ +TEST_CASE( "AO-system configuration reports atmosphere validity", "[ao::analysis::aoSystem]" ) +{ + aoSystem> aoSystem; + appConfigurator config; + aoSystem.setupConfig( config ); + + SECTION( "default atmosphere" ) + { + REQUIRE( aoSystem.loadConfig( config ) == mx::error_t::sizeerr ); + } + + SECTION( "complete atmosphere" ) + { + aoSystem.atm.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + REQUIRE( aoSystem.loadConfig( config ) == mx::error_t::noerror ); + } +} 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/ao/analysis/clGainOpt_test.cpp b/tests/include/ao/analysis/clGainOpt_test.cpp new file mode 100644 index 000000000..dccc01ba3 --- /dev/null +++ b/tests/include/ao/analysis/clGainOpt_test.cpp @@ -0,0 +1,224 @@ +/** \file clGainOpt_test.cpp + * \brief Tests of closed-loop gain optimization transfer functions. + */ + +#include "../../../catch2/catch.hpp" + +#define MX_NO_ERROR_REPORTS + +#include "../../../../include/ao/analysis/clGainOpt.hpp" + +#include +#include +#include +#include + +namespace +{ + +using optimizerT = mx::AO::analysis::clGainOpt; + +/// Require two complex values to agree within floating-point precision. +void requireComplexEqual( const std::complex &actual, /**< [in] value produced by the retimed optimizer */ + const std::complex &expected /**< [in] value produced by the fresh optimizer */ ) +{ + REQUIRE( actual.real() == Approx( expected.real() ).epsilon( 1e-12 ).margin( 1e-14 ) ); + REQUIRE( actual.imag() == Approx( expected.imag() ).epsilon( 1e-12 ).margin( 1e-14 ) ); +} + +} // namespace + +/// Verify that changing the sampling interval invalidates all sampling-interval-dependent cached values. +/** Exercises mx::AO::analysis::clGainOpt::Ti and the public transfer-function calculations after the trigonometric + * cache has already been populated. + */ +TEST_CASE( "Gain optimizer recomputes transfer functions after changing Ti", "[ao::analysis::clGainOpt]" ) +{ + constexpr double initialTi = 0.001; + constexpr double updatedTi = 0.0017; + constexpr double delay = 0.0013; + constexpr double gain = 0.42; + + const std::vector frequency{ 25.0, 50.0, 75.0, 100.0 }; + const std::vector fir{ 0.75, 0.20, -0.05 }; + const std::vector iir{ 0.45, -0.08, 0.025 }; + + optimizerT retimed( initialTi, delay ); + retimed.f( frequency ); + retimed.b( fir ); + retimed.a( iir ); + retimed.remember( 0.97 ); + + for( std::size_t index = 0; index < frequency.size(); ++index ) + { + retimed.olXfer( index ); + } + + retimed.Ti( updatedTi ); + + optimizerT fresh( updatedTi, delay ); + fresh.f( frequency ); + fresh.b( fir ); + fresh.a( iir ); + fresh.remember( 0.97 ); + + for( std::size_t index = 0; index < frequency.size(); ++index ) + { + CAPTURE( index, frequency[index] ); + + optimizerT::complexT retimedDm; + optimizerT::complexT retimedDelay; + optimizerT::complexT retimedController; + const optimizerT::complexT retimedOpenLoop = + retimed.olXfer( index, retimedDm, retimedDelay, retimedController ); + + optimizerT::complexT freshDm; + optimizerT::complexT freshDelay; + optimizerT::complexT freshController; + const optimizerT::complexT freshOpenLoop = fresh.olXfer( index, freshDm, freshDelay, freshController ); + + requireComplexEqual( retimedDm, freshDm ); + requireComplexEqual( retimedDelay, freshDelay ); + requireComplexEqual( retimedController, freshController ); + requireComplexEqual( retimedOpenLoop, freshOpenLoop ); + requireComplexEqual( retimed.clETF( index, gain ), fresh.clETF( index, gain ) ); + REQUIRE( retimed.clETFPhase( index, gain ) == + Approx( fresh.clETFPhase( index, gain ) ).epsilon( 1e-12 ).margin( 1e-14 ) ); + REQUIRE( retimed.clETF2( index, gain ) == + Approx( fresh.clETF2( index, gain ) ).epsilon( 1e-12 ).margin( 1e-14 ) ); + requireComplexEqual( retimed.clNTF( index, gain ), fresh.clNTF( index, gain ) ); + REQUIRE( retimed.clNTF2( index, gain ) == + Approx( fresh.clNTF2( index, gain ) ).epsilon( 1e-12 ).margin( 1e-14 ) ); + + double retimedEtf = 0; + double retimedNtf = 0; + retimed.clTF2( retimedEtf, retimedNtf, index, gain ); + + double freshEtf = 0; + double freshNtf = 0; + fresh.clTF2( freshEtf, freshNtf, index, gain ); + + REQUIRE( retimedEtf == Approx( freshEtf ).epsilon( 1e-12 ).margin( 1e-14 ) ); + REQUIRE( retimedNtf == Approx( freshNtf ).epsilon( 1e-12 ).margin( 1e-14 ) ); + } + + const std::vector disturbancePsd{ 4.0, 2.0, 0.7, 0.2 }; + const std::vector noisePsd{ 0.01, 0.01, 0.01, 0.01 }; + REQUIRE( retimed.clVariance( disturbancePsd, noisePsd, gain ) == + Approx( fresh.clVariance( disturbancePsd, noisePsd, gain ) ).epsilon( 1e-12 ).margin( 1e-14 ) ); + + double retimedOptimalVariance = 0; + double retimedOptimalGain = 0; + REQUIRE( + retimed.optGainOpenLoop( retimedOptimalGain, retimedOptimalVariance, disturbancePsd, noisePsd, 0.5, false ) == + mx::error_t::noerror ); + + double freshOptimalVariance = 0; + double freshOptimalGain = 0; + REQUIRE( fresh.optGainOpenLoop( freshOptimalGain, freshOptimalVariance, disturbancePsd, noisePsd, 0.5, false ) == + mx::error_t::noerror ); + + REQUIRE( retimedOptimalGain == Approx( freshOptimalGain ).epsilon( 1e-12 ).margin( 1e-14 ) ); + REQUIRE( retimedOptimalVariance == Approx( freshOptimalVariance ).epsilon( 1e-12 ).margin( 1e-14 ) ); +} + +/// Verify that maximum-stable-gain calculation interpolates a bracketed pure-integrator Nyquist crossing. +/** Exercises mx::AO::analysis::clGainOpt::maxStableGain and its crossing diagnostics. */ +TEST_CASE( "Gain optimizer interpolates a pure-integrator stability crossing", "[ao::analysis::clGainOpt]" ) +{ + optimizerT optimizer( 0.001, 0.0015 ); + const std::vector frequency{ 100.0, 124.0, 126.0, 150.0 }; + optimizer.f( frequency ); + + const double crossingPhase = std::numbers::pi_v / 4.0; + const double expectedGain = crossingPhase * crossingPhase / ( 2.0 * std::sin( crossingPhase / 2.0 ) ); + const double lowerSampleGain = -1.0 / optimizer.olXfer( 1 ).real(); + + double maximumGain = 0; + optimizerT::maxStableGainReport report; + REQUIRE( optimizer.maxStableGain( maximumGain, &report ) == mx::error_t::noerror ); + + REQUIRE( report.status == optimizerT::maxStableGainStatus::crossingFound ); + REQUIRE( report.lowerIndex == 1 ); + REQUIRE( report.upperIndex == 2 ); + REQUIRE( report.lowerFrequency == 124.0 ); + REQUIRE( report.upperFrequency == 126.0 ); + REQUIRE( report.crossingFrequency == Approx( 125.0 ).margin( 0.02 ) ); + REQUIRE( maximumGain == report.gain ); + REQUIRE( std::abs( maximumGain - expectedGain ) < std::abs( lowerSampleGain - expectedGain ) ); +} + +/// Verify that maximum-stable-gain calculation reports invalid grids and missing crossings. +/** Exercises mx::AO::analysis::clGainOpt::maxStableGain failure statuses without sentinel gain values. */ +TEST_CASE( "Gain optimizer reports missing stability crossings", "[ao::analysis::clGainOpt]" ) +{ + optimizerT optimizer( 0.001, 0.0015 ); + double maximumGain = 0; + optimizerT::maxStableGainReport report; + + optimizer.f( std::vector{ 1.0 } ); + REQUIRE( optimizer.maxStableGain( maximumGain, &report ) == mx::error_t::sizeerr ); + REQUIRE( report.status == optimizerT::maxStableGainStatus::invalidInput ); + REQUIRE( mx::math::isNan( maximumGain ) ); + + optimizer.f( std::vector{ 1.0, 2.0, 3.0 } ); + REQUIRE( optimizer.maxStableGain( maximumGain, &report ) == mx::error_t::notfound ); + REQUIRE( report.status == optimizerT::maxStableGainStatus::noCrossing ); + REQUIRE( mx::math::isNan( maximumGain ) ); +} + +/// Verify that optimum-gain calculation enforces its interval and reports termination state. +/** Exercises both mx::AO::analysis::clGainOpt::optGainOpenLoop overloads for small intervals, invalid intervals, + * stability failure, and forced iteration exhaustion. + */ +TEST_CASE( "Gain optimizer reports minimizer termination", "[ao::analysis::clGainOpt]" ) +{ + optimizerT optimizer( 0.001, 0.0015 ); + std::vector frequency; + std::vector disturbance; + std::vector noise; + for( int index = 1; index <= 500; ++index ) + { + const double value = static_cast( index ); + frequency.push_back( value ); + disturbance.push_back( 1.0 / ( 1.0 + value * value ) ); + noise.push_back( 1e-6 ); + } + optimizer.f( frequency ); + + double optimalGain = 0; + double variance = 0; + optimizerT::optGainReport report; + constexpr double smallMaximumGain = 1e-3; + REQUIRE( optimizer.optGainOpenLoop( optimalGain, variance, disturbance, noise, smallMaximumGain, true, &report ) == + mx::error_t::noerror ); + REQUIRE( ( report.status == optimizerT::optGainStatus::converged || + report.status == optimizerT::optGainStatus::boundaryLimited ) ); + REQUIRE( report.minimumEvaluatedGain >= optimizer.m_minFindMin ); + REQUIRE( report.maximumEvaluatedGain <= optimizer.m_minFindMaxFact * smallMaximumGain ); + REQUIRE( mx::math::isFinite( optimalGain ) ); + REQUIRE( mx::math::isFinite( variance ) ); + + REQUIRE( optimizer.optGainOpenLoop( optimalGain, variance, disturbance, noise, 1e-10, false, &report ) == + mx::error_t::invalidconfig ); + REQUIRE( report.status == optimizerT::optGainStatus::invalidInput ); + REQUIRE( mx::math::isNan( optimalGain ) ); + REQUIRE( mx::math::isNan( variance ) ); + + optimizer.f( std::vector{ 1.0, 2.0, 3.0 } ); + disturbance.resize( 3 ); + noise.resize( 3 ); + REQUIRE( optimizer.optGainOpenLoop( optimalGain, variance, disturbance, noise, false, &report ) == + mx::error_t::notfound ); + REQUIRE( report.status == optimizerT::optGainStatus::stabilityFailure ); + REQUIRE( report.stability.status == optimizerT::maxStableGainStatus::noCrossing ); + + optimizer.f( frequency ); + disturbance.resize( frequency.size(), 1e-3 ); + noise.resize( frequency.size(), 1e-6 ); + optimizer.m_minFindMaxIter = 1; + REQUIRE( optimizer.optGainOpenLoop( optimalGain, variance, disturbance, noise, 0.5, false, &report ) == + mx::error_t::timeout ); + REQUIRE( report.status == optimizerT::optGainStatus::iterationLimit ); + REQUIRE( report.iterations == 1 ); +} diff --git a/tests/include/ao/analysis/fourierTemporalPSD_test.cpp b/tests/include/ao/analysis/fourierTemporalPSD_test.cpp new file mode 100644 index 000000000..2ae35ae7a --- /dev/null +++ b/tests/include/ao/analysis/fourierTemporalPSD_test.cpp @@ -0,0 +1,555 @@ +/** \file fourierTemporalPSD_test.cpp + * \brief Tests of Fourier-mode temporal power spectral densities. + */ + +#include "../../../catch2/catch.hpp" + +#define MX_NO_ERROR_REPORTS + +#include "../../../../include/ao/analysis/fourierTemporalPSD.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using aoSystemBaseT = mx::AO::analysis::aoSystem>; + +/// AO-system test type that forces local template instantiation under sanitizers. +struct aoSystemT : public aoSystemBaseT +{ +}; + +using temporalPsdT = mx::AO::analysis::fourierTemporalPSD; + +namespace +{ +/// \cond fourierTemporalPSD_test_detail + +/// Sentinel GSL error handler used to verify restoration after PSD calculations. +void testGslErrorHandler( const char *reason, /**< [in] GSL diagnostic text */ + const char *file, /**< [in] GSL source file */ + int line, /**< [in] GSL source line */ + int errorNumber /**< [in] GSL status code */ ) +{ + static_cast( reason ); + static_cast( file ); + static_cast( line ); + static_cast( errorNumber ); +} + +/// Return a null workspace to exercise allocation-failure propagation. +gsl_integration_workspace *failWorkspaceAllocation( size_t size /**< [in] requested workspace size */ ) +{ + static_cast( size ); + return nullptr; +} + +/// Test evaluator exposing the protected allocator-injection constructor. +struct allocationFailureTemporalPsdT : public temporalPsdT +{ + /// Construct an evaluator whose workspace allocation always fails. + allocationFailureTemporalPsdT() : temporalPsdT( &failWorkspaceAllocation ) + { + } +}; + +/// Test evaluator exposing workspace allocation and ownership state. +struct workspaceOwnershipTemporalPsdT : public temporalPsdT +{ + /// Preserve move construction for the derived test evaluator. + workspaceOwnershipTemporalPsdT( workspaceOwnershipTemporalPsdT && ) noexcept = default; + + /// Construct a test evaluator without allocating a workspace. + workspaceOwnershipTemporalPsdT() = default; + + /// Allocate the workspace through the production helper. + mx::error_t allocateTestWorkspace() + { + return allocateWorkspace(); + } + + /// Return whether this evaluator currently owns a workspace. + bool ownsWorkspace() const + { + return m_workspace != nullptr; + } +}; + +/// \endcond + +} // namespace + +/// Verify that Fourier temporal PSD evaluators uniquely own movable GSL workspaces. +/** Exercises the ownership contract of mx::AO::analysis::fourierTemporalPSD. */ +TEST_CASE( "Fourier temporal PSD workspace ownership is move-only", "[ao::analysis::fourierTemporalPSD]" ) +{ + STATIC_REQUIRE_FALSE( std::is_copy_constructible_v ); + STATIC_REQUIRE_FALSE( std::is_copy_assignable_v ); + STATIC_REQUIRE( std::is_nothrow_move_constructible_v ); + STATIC_REQUIRE( std::is_nothrow_move_assignable_v ); + + workspaceOwnershipTemporalPsdT source; + REQUIRE( source.allocateTestWorkspace() == mx::error_t::noerror ); + REQUIRE( source.ownsWorkspace() ); + + workspaceOwnershipTemporalPsdT destination( std::move( source ) ); + REQUIRE_FALSE( source.ownsWorkspace() ); + REQUIRE( destination.ownsWorkspace() ); + REQUIRE( destination.m_aosys == nullptr ); +} + +/// Verify that public Fourier temporal PSD calculations reject malformed inputs before modifying output. +/** Exercises mx::AO::analysis::fourierTemporalPSD::singleLayerPSD and + * mx::AO::analysis::fourierTemporalPSD::multiLayerPSD precondition handling. */ +TEST_CASE( "Fourier temporal PSD validates public calculation inputs", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.lam_sci( 1.0e-6 ); + aoSystem.lam_wfs( 0.8e-6 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + std::vector frequency{ 0.0, 1.0 }; + std::vector psd( frequency.size(), 7.0 ); + const std::vector unchanged = psd; + temporalPsdT::reportT report; + report.record( GSL_SUCCESS, 0, 0.0, 0.0, 0.0, 1.0, 1.0 ); + + SECTION( "null AO system" ) + { + temporalPsd.m_aosys = nullptr; + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, frequency.back(), &report ) == + mx::error_t::invalidconfig ); + } + + SECTION( "empty vectors" ) + { + std::vector empty; + REQUIRE( temporalPsd.singleLayerPSD( empty, empty, 1.0, 0.0, 0, 1, 0, &report ) == mx::error_t::sizeerr ); + } + + SECTION( "mismatched vector sizes" ) + { + std::vector undersizedPsd{ 7.0 }; + REQUIRE( temporalPsd.singleLayerPSD( undersizedPsd, frequency, 1.0, 0.0, 0, 1, frequency.back(), &report ) == + mx::error_t::sizeerr ); + REQUIRE( undersizedPsd == std::vector{ 7.0 } ); + } + + SECTION( "nonmonotone frequency grid" ) + { + frequency[1] = frequency[0]; + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == mx::error_t::invalidarg ); + } + + SECTION( "negative frequency" ) + { + frequency[0] = -1.0; + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == mx::error_t::invalidarg ); + } + + SECTION( "nonfinite frequency grid" ) + { + frequency[1] = std::numeric_limits::quiet_NaN(); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == mx::error_t::invalidarg ); + } + + SECTION( "nonfinite mode" ) + { + REQUIRE( + temporalPsd + .singleLayerPSD( psd, frequency, std::numeric_limits::infinity(), 0.0, 0, 1, 0, &report ) == + mx::error_t::invalidarg ); + } + + SECTION( "invalid parity" ) + { + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 0, 0, &report ) == mx::error_t::invalidarg ); + } + + SECTION( "negative cutoff" ) + { + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, -1.0, &report ) == + mx::error_t::invalidarg ); + } + + SECTION( "invalid layer index" ) + { + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 1, 1, 0, &report ) == mx::error_t::invalidarg ); + } + + SECTION( "invalid tolerance" ) + { + temporalPsd.absTol( 0 ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == + mx::error_t::invalidconfig ); + } + + SECTION( "invalid relative tolerance" ) + { + temporalPsd.relTol( 1 ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == + mx::error_t::invalidconfig ); + } + + SECTION( "nonpositive aperture" ) + { + aoSystem.D( 0 ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == + mx::error_t::invalidconfig ); + } + + SECTION( "zero wind" ) + { + aoSystem.atm.layer_v_wind( std::vector{ 0.0 } ); + REQUIRE( temporalPsd.multiLayerPSD( psd, frequency, 1.0, 0.0, 1, 0, &report ) == + mx::error_t::invalidconfig ); + } + + SECTION( "mismatched atmosphere vectors" ) + { + aoSystem.atm.layer_dir( std::vector{} ); + REQUIRE( temporalPsd.multiLayerPSD( psd, frequency, 1.0, 0.0, 1, 0, &report ) == mx::error_t::sizeerr ); + } + + REQUIRE( psd == unchanged ); + REQUIRE( report.integrationsAttempted == 0 ); +} + +/// Verify that GSL workspace allocation failure is returned without evaluating or modifying the PSD. +/** Exercises mx::AO::analysis::fourierTemporalPSD::singleLayerPSD with an injected failing workspace allocator. */ +TEST_CASE( "Fourier temporal PSD reports workspace allocation failure", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + allocationFailureTemporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + std::vector frequency{ 1.0 }; + std::vector psd{ 7.0 }; + temporalPsdT::reportT report; + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, 0, &report ) == mx::error_t::allocerr ); + REQUIRE( psd == std::vector{ 7.0 } ); + REQUIRE( report.integrationsAttempted == 0 ); +} + +/// Verify that PSD-grid generation validates its public inputs before creating output. +/** Exercises mx::AO::analysis::fourierTemporalPSD::makePSDGrid precondition handling. */ +TEST_CASE( "Fourier temporal PSD validates grid-generation inputs", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + const std::filesystem::path outputDirectory = + std::filesystem::temp_directory_path() / + ( "mxlib_fourierTemporalPSD_invalid_grid_" + std::to_string( static_cast( getpid() ) ) ); + std::error_code filesystemError; + std::filesystem::remove_all( outputDirectory, filesystemError ); + REQUIRE_FALSE( filesystemError ); + + SECTION( "empty output directory" ) + { + REQUIRE( temporalPsd.makePSDGrid( "", 1, 1.0, 1.0 ) == mx::error_t::invalidarg ); + } + + SECTION( "nonpositive spatial extent" ) + { + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 0, 1.0, 1.0 ) == mx::error_t::invalidarg ); + } + + SECTION( "nonpositive frequency spacing" ) + { + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 1, 0.0, 1.0 ) == mx::error_t::invalidarg ); + } + + SECTION( "nonfinite frequency spacing" ) + { + REQUIRE( + temporalPsd.makePSDGrid( outputDirectory.string(), 1, std::numeric_limits::quiet_NaN(), 1.0 ) == + mx::error_t::invalidarg ); + } + + SECTION( "nonpositive maximum frequency" ) + { + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 1, 1.0, 0.0 ) == mx::error_t::invalidarg ); + } + + SECTION( "negative exact-calculation cutoff" ) + { + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 1, 1.0, 1.0, -1.0 ) == mx::error_t::invalidarg ); + } + + SECTION( "unrepresentable sample count" ) + { + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 1, std::numeric_limits::min(), 1.0 ) == + mx::error_t::sizeerr ); + } + + SECTION( "invalid AO system" ) + { + temporalPsd.m_aosys = nullptr; + REQUIRE( temporalPsd.makePSDGrid( outputDirectory.string(), 1, 1.0, 1.0 ) == mx::error_t::invalidconfig ); + } + + REQUIRE_FALSE( std::filesystem::exists( outputDirectory ) ); +} + +/// Verify that PSD-grid generation propagates calculation and output failures. +/** Exercises mx::AO::analysis::fourierTemporalPSD::makePSDGrid failure handling after successful preflight. */ +TEST_CASE( "Fourier temporal PSD reports grid-generation failures", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + const std::filesystem::path outputPath = + std::filesystem::temp_directory_path() / + ( "mxlib_fourierTemporalPSD_grid_failure_" + std::to_string( static_cast( getpid() ) ) ); + std::error_code filesystemError; + std::filesystem::remove_all( outputPath, filesystemError ); + REQUIRE_FALSE( filesystemError ); + + SECTION( "mode calculation failure" ) + { + allocationFailureTemporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + REQUIRE( temporalPsd.makePSDGrid( outputPath.string(), 1, 1.0, 1.0 ) == mx::error_t::allocerr ); + REQUIRE( std::filesystem::exists( outputPath / "params.txt" ) ); + REQUIRE( std::filesystem::exists( outputPath / "psds" / "freq.binv" ) ); + } + + SECTION( "output path is a regular file" ) + { + std::ofstream outputFile( outputPath ); + REQUIRE( outputFile.is_open() ); + outputFile.close(); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + REQUIRE( temporalPsd.makePSDGrid( outputPath.string(), 1, 1.0, 1.0 ) == mx::error_t::enotdir ); + } + + std::filesystem::remove_all( outputPath, filesystemError ); + REQUIRE_FALSE( filesystemError ); +} + +/// Verify that Fourier PSD calculation requires a valid atmosphere and accepts every complete preset. +/** Exercises mx::AO::analysis::aoAtmosphere::validate, + * mx::AO::analysis::aoAtmosphere::loadGuyon2005, mx::AO::analysis::aoAtmosphere::loadLCO, + * mx::AO::analysis::aoAtmosphere::setSingleLayer, and + * mx::AO::analysis::fourierTemporalPSD::singleLayerPSD. */ +TEST_CASE( "Fourier temporal PSD enforces atmosphere validity at calculation", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + std::vector frequency{ 1.0 }; + std::vector psd{ 7.0 }; + + SECTION( "default atmosphere" ) + { + REQUIRE( aoSystem.atm.validate() == mx::error_t::sizeerr ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, frequency.back() ) == + mx::error_t::sizeerr ); + REQUIRE( psd == std::vector{ 7.0 } ); + } + + SECTION( "Guyon 2005 preset with infinite outer scale" ) + { + aoSystem.atm.loadGuyon2005(); + REQUIRE( aoSystem.atm.L_0( 0 ) == 0.0 ); + REQUIRE( aoSystem.atm.validate() == mx::error_t::noerror ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, frequency.back() ) == + mx::error_t::noerror ); + } + + SECTION( "LCO preset" ) + { + aoSystem.atm.loadLCO(); + REQUIRE( aoSystem.atm.validate() == mx::error_t::noerror ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, frequency.back() ) == + mx::error_t::noerror ); + } + + SECTION( "single-layer preset" ) + { + aoSystem.atm.setSingleLayer( 0.2, 0.5e-6, 25.0, 0.0, 0.0, 10.0, 0.0 ); + REQUIRE( aoSystem.atm.validate() == mx::error_t::noerror ); + REQUIRE( temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, frequency.back() ) == + mx::error_t::noerror ); + } +} + +/// Verify temporal-PSD tail initialization at and below its nominal averaging width. +/** Exercises zero, one, 49, and 50 exactly integrated frequency bins. */ +TEST_CASE( "Fourier temporal PSD handles short exact tails", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = mx::AO::analysis::basis::basic; + + constexpr std::size_t frequencyCount = 55; + std::vector frequency( frequencyCount ); + for( std::size_t index = 0; index < frequency.size(); ++index ) + { + frequency[index] = static_cast( index + 1 ); + } + + for( const std::size_t exactCount : { 0U, 1U, 49U, 50U } ) + { + DYNAMIC_SECTION( exactCount << " exact bins" ) + { + std::vector psd( frequency.size(), 0.0 ); + const double maximumExactFrequency = exactCount == 0 ? 0.5 : frequency[exactCount - 1]; + temporalPsdT::reportT report; + + const mx::error_t result = + temporalPsd.singleLayerPSD( psd, frequency, 1.0, 0.0, 0, 1, maximumExactFrequency, &report ); + + if( exactCount == 0 ) + { + REQUIRE( result == mx::error_t::invalidarg ); + REQUIRE( report.integrationsAttempted == 0 ); + continue; + } + + REQUIRE( result == mx::error_t::noerror ); + REQUIRE( report.integrationsAttempted == exactCount ); + + const std::size_t averageCount = std::min( exactCount, 50 ); + const double exponent = aoSystem.atm.alpha( 0 ) + 2.0; + double expected = 0.0; + for( std::size_t offset = averageCount; offset > 0; --offset ) + { + const std::size_t index = exactCount - offset; + expected += psd[index] * std::pow( frequency[index] / frequency[exactCount], exponent ); + } + expected /= static_cast( averageCount ); + + REQUIRE( psd[exactCount] == Approx( expected ).epsilon( 1e-12 ) ); + REQUIRE( std::isfinite( psd[exactCount] ) ); + REQUIRE( psd[exactCount] >= 0.0 ); + + const double expectedLast = + psd[exactCount] * std::pow( frequency[exactCount] / frequency.back(), exponent ); + REQUIRE( psd.back() == Approx( expectedLast ).epsilon( 1e-12 ) ); + } + } +} + +/// Verify aggregation and formatting of Fourier temporal-PSD quadrature diagnostics. +/** Exercises `fourierTemporalPSDReport::record`, `fourierTemporalPSDReport::merge`, and + * `fourierTemporalPSDReport::write` with multiple GSL statuses and atmospheric layers. */ +TEST_CASE( "Fourier temporal PSD aggregates quadrature diagnostics", "[ao::analysis::fourierTemporalPSD]" ) +{ + temporalPsdT::reportT report; + report.record( GSL_SUCCESS, 0, 1.0, 2.0, 0.1, 0.1, 0.01 ); + report.record( GSL_EROUND, 0, 2.0, 2.0, 0.4, 0.1, 0.01 ); + report.record( GSL_EROUND, 1, 3.0, 4.0, 1.0, 0.1, 0.01 ); + + temporalPsdT::reportT other; + other.record( GSL_EDIVERGE, 2, 4.0, 1.0, 0.5, 0.1, 0.01 ); + report.merge( other ); + + REQUIRE( report.integrationsAttempted == 4 ); + REQUIRE( report.integrationsConverged == 1 ); + REQUIRE( report.failureCount() == 3 ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).count == 2 ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).countByLayer.at( 0 ) == 1 ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).countByLayer.at( 1 ) == 1 ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).maximumAbsoluteError == Approx( 1.0 ) ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).maximumToleranceRatio == Approx( 10.0 ) ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).worstLayer == 1 ); + REQUIRE( report.gslStatus.at( GSL_EROUND ).worstFrequency == Approx( 3.0 ) ); + REQUIRE( report.gslStatus.at( GSL_EDIVERGE ).count == 1 ); + + std::ostringstream summary; + report.write( summary ); + REQUIRE( summary.str().find( "1/4 converged" ) != std::string::npos ); + REQUIRE( summary.str().find( gsl_strerror( GSL_EROUND ) ) != std::string::npos ); + REQUIRE( summary.str().find( gsl_strerror( GSL_EDIVERGE ) ) != std::string::npos ); + REQUIRE( summary.str().find( "layers {0: 1, 1: 1}" ) != std::string::npos ); +} + +/// Verify permissive and strict handling of GSL quadrature statuses. +/** Exercises the status policy used by `fourierTemporalPSD::singleLayerPSD`. */ +TEST_CASE( "Fourier temporal PSD applies quadrature policy", "[ao::analysis::fourierTemporalPSD]" ) +{ + using mx::AO::analysis::fourierTemporalPSDPolicy; + using mx::AO::analysis::fourierTemporalPSD_detail::applyPolicy; + + REQUIRE( applyPolicy( GSL_SUCCESS, fourierTemporalPSDPolicy::permissive ) == mx::error_t::noerror ); + for( const int status : { GSL_EMAXITER, GSL_EROUND, GSL_ESING, GSL_EDIVERGE } ) + { + REQUIRE( applyPolicy( status, fourierTemporalPSDPolicy::permissive ) == mx::error_t::noerror ); + REQUIRE( applyPolicy( status, fourierTemporalPSDPolicy::strict ) == mx::error_t::liberr ); + } + REQUIRE( applyPolicy( GSL_EDOM, fourierTemporalPSDPolicy::permissive ) == mx::error_t::invalidconfig ); + REQUIRE( applyPolicy( GSL_EINVAL, fourierTemporalPSDPolicy::permissive ) == mx::error_t::invalidconfig ); + REQUIRE( applyPolicy( GSL_ENOMEM, fourierTemporalPSDPolicy::permissive ) == mx::error_t::allocerr ); + REQUIRE( applyPolicy( GSL_EFAILED, fourierTemporalPSDPolicy::permissive ) == mx::error_t::liberr ); +} + +/// Verify multilayer error propagation and restoration of the caller's GSL handler. +/** Exercises `fourierTemporalPSD::multiLayerPSD` with an invalid basis. */ +TEST_CASE( "Fourier temporal PSD propagates multilayer errors", "[ao::analysis::fourierTemporalPSD]" ) +{ + aoSystemT aoSystem; + aoSystem.D( 6.5 ); + aoSystem.atm.setSingleLayer( 0.16, 500e-9, 25.0, 0.0, 0.0, 10.0, 0.0 ); + + temporalPsdT temporalPsd; + temporalPsd.m_aosys = &aoSystem; + temporalPsd._useBasis = -1; + + std::vector frequency{ 1.0 }; + std::vector psd( frequency.size(), 0.0 ); + temporalPsdT::reportT report; + + gsl_error_handler_t *previousHandler = gsl_set_error_handler( &testGslErrorHandler ); + const mx::error_t result = + temporalPsd.multiLayerPSD( psd, + frequency, + 1.0, + 0.0, + 1, + frequency.back(), + &report, + mx::AO::analysis::fourierTemporalPSDPolicy::permissive ); + gsl_error_handler_t *observedHandler = gsl_set_error_handler( previousHandler ); + + REQUIRE( result == mx::error_t::invalidarg ); + REQUIRE( report.integrationsAttempted == 0 ); + REQUIRE( observedHandler == &testGslErrorHandler ); +} 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() ) ); +}