clAOLinearPredictor.hpp

Go to the documentation of this file.

/** \file clAOLinearPredictor.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Provides a class to manage closed loop gain linear predictor determination.
 * \ingroup mxAO_files
 *
 */
 
#ifndef clAOLinearPredictor_hpp
#define clAOLinearPredictor_hpp
 
#include <vector>
 
#include "../../math/geo.hpp"
 
#include "../../sigproc/psdUtils.hpp"
 
#include "../../sigproc/autocorrelation.hpp"
#include "../../sigproc/linearPredictor.hpp"
 
#include "clGainOpt.hpp"
 
namespace mx
{
namespace AO
{
namespace analysis
{
 
/// Class to manage the calculation of linear predictor coefficients for a closed-loop AO system.
/**
 * \tparam _realT the real floating point type in which to do all arithmetic.
 *
 * \ingroup mxAOAnalytic
 */
template <typename _realT>
struct clAOLinearPredictor
{
    typedef _realT realT;
 
    std::vector<realT> PSDtn;
 
    std::vector<realT> psd2s;
 
    std::vector<realT> ac;
 
    sigproc::autocorrelationFromPSD<realT> acpsd;
 
    sigproc::linearPredictor<realT> lp;
 
    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_gmax_lp{ 5 }; ///< The maximum allowable gain for LP.
 
    // Stopping conditions:
    realT m_minPrecision{ 0.001 };
    int m_maxIts{ 100 };
 
    int extrap;
 
    clAOLinearPredictor()
    {
        extrap = 0;
    }
 
    /// 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.
     *
     */
    int calcCoefficients( std::vector<realT> &PSDt, ///< [in] the turbulence PSD
                          std::vector<realT> &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 )
    {
        PSDtn.resize( PSDt.size() );
 
        for( size_t i = 0; i < PSDt.size(); ++i )
        {
            PSDtn[i] = PSDt[i] + PSDn[i] + PSDreg;
        }
 
        sigproc::augment1SidedPSD( psd2s, PSDtn, 1 );
 
        ac.resize( psd2s.size() );
 
        acpsd( ac, psd2s );
 
        return lp.calcCoefficients( ac, Nc, condition, extrap );
    }
 
    /// Worker function for regularizing the PSD for coefficient calculation.
    /**
     * \tparam printout if true then the results are printed to stdout as they are calculated.
     *
     * On first call (min_var = 0):
     *     loop over scale factors from min_sc to max_sc (<=) in steps of precision.
     *
     * 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
     */
    template <bool printout>
    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<realT> &go_lp,  ///< [in] the gain optimization object
                                 std::vector<realT> &PSDt, ///< [in] the turbulence PSD
                                 std::vector<realT> &PSDn, ///< [in] the WFS noise PSD
                                 int Nc                    ///< [in] the number of coefficients
    )
    {
        realT gmax_lp;
        realT gopt_lp;
        realT var_lp;
 
        realT sc0;
 
        // realT last_var;
 
        // min_var == 0 indicates first call
        if( min_var == 0 )
        {
            sc0 = min_sc;
            min_var = std::numeric_limits<realT>::max();
            // last_var = std::numeric_limits<realT>::max();
        }
        else
        {
            sc0 = min_sc - precision * m_dPrecision;
            // sc0 = min_sc;
            // last_var = min_var;
        }
 
        // auto it = std::max_element(std::begin(PSDt), std::end(PSDt));
        realT psdReg = PSDt[0]; //*it/10;
 
        // Test from sc0 to max_sc in steps of precision
        for( realT sc = sc0; sc <= max_sc; sc += precision )
        {
            if( calcCoefficients( PSDt, PSDn, psdReg * pow( 10, -sc / 10 ), Nc ) < 0 )
                return -1;
 
            go_lp.a( lp.m_c );
            go_lp.b( lp.m_c );
 
            realT ll = 0, ul = 0;
            gmax_lp = go_lp.maxStableGain( ll, ul );
            if( gmax_lp > m_gmax_lp )
            {
                gmax_lp = m_gmax_lp;
            }
 
            gopt_lp = go_lp.optGainOpenLoop( var_lp, PSDt, PSDn, gmax_lp, false );
 
            if( printout )
            {
                std::cout << -( sc ) / 10 << " " << gmax_lp << " " << gopt_lp << " " << var_lp << "\n";
            }
 
            if( var_lp < min_var )
            {
                min_var = var_lp;
                min_sc = sc;
            }
 
            // A jump by a factor of 10 indicates the wall
            if( var_lp > 10 * min_var )
                return 0;
 
            // last_var = var_lp;
        }
 
        return -1;
    }
 
    /// 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.
     */
    template <bool printout = false>
    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.
                                clGainOpt<realT> &go_lp,  ///< [in] the gain optimization object
                                std::vector<realT> &PSDt, ///< [in] the turbulence PSD
                                std::vector<realT> &PSDn, ///< [in] the WFS noise PSD
                                int Nc                    ///< [in] the number of coefficients
    )
    {
 
        realT min_var = m_min_var0;
        realT min_sc = m_min_sc0;
        realT precision = m_precision0;
        realT max_sc = m_max_sc0;
 
        int its = 0;
        while( precision > m_minPrecision && its < m_maxIts )
        {
            _regularizeCoefficients<printout>( min_var, min_sc, precision, max_sc, go_lp, PSDt, PSDn, Nc );
 
            if( min_sc == max_sc )
            {
                if( its == 0 )
                {
                    min_sc -= precision;
                    max_sc = 200;
                }
                else
                {
                    // std::cerr << "Error in regularizeCoefficients.\n";
                    // return -1;
                    break;
                }
            }
            else
            {
                max_sc = min_sc + precision;
                precision /= m_dPrecision;
            }
 
            ++its;
        }
 
        // Now record final values
        if( calcCoefficients( PSDt, PSDn, PSDt[0] * pow( 10, -min_sc / 10 ), Nc ) < 0 )
            return -1;
 
        go_lp.a( lp.m_c );
        // go_lp.a(std::vector<realT>({1}));
        go_lp.b( lp.m_c );
        // go_lp.b(std::vector<realT>({lp._c(0,0)}));
 
        realT ll = 0, ul = 0;
        gmax_lp = go_lp.maxStableGain( ll, ul );
        gopt_lp = go_lp.optGainOpenLoop( var_lp, PSDt, PSDn, gmax_lp, false );
 
        return 0;
    }
 
    /// 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.
     */
    template <bool printout = false>
    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<realT> &go_lp,  ///< [in] the gain optimization object
                    std::vector<realT> &PSDt, ///< [in] the turbulence PSD
                    std::vector<realT> &PSDn, ///< [in] the WFS noise PSD
                    int minNc,                ///< [in] the number of coefficients
                    int maxNc )
    {
        realT minVar = std::numeric_limits<realT>::max();
 
        for( int n = minNc; n <= maxNc; ++n )
        {
            realT _gmax_lp;
            realT _gopt_lp;
            realT _var_lp;
            regularizeCoefficients<printout>( _gmax_lp, _gopt_lp, _var_lp, go_lp, PSDt, PSDn, n );
 
            if( _var_lp < minVar )
            {
                gmax_lp = _gmax_lp;
                gopt_lp = _gopt_lp;
                var_lp = _var_lp;
                Nc = n;
 
                minVar = var_lp;
            }
        }
 
        return 0;
    }
};
 
} // namespace analysis
} // namespace AO
} // namespace mx
 
#endif // clAOLinearPredictor_hpp