Vectorfit.h

#ifndef CVectorfit_h__
#define CVectorfit_h__
 
#include <complex>
#include <Eigen/Eigen>
 
namespace Ship{
    typedef Eigen::Matrix< std::complex<double>,1,Eigen::Dynamic> vec_complex;
    typedef Eigen::Matrix< double,1,Eigen::Dynamic> vec;
    typedef Eigen::Matrix< double,Eigen::Dynamic,Eigen::Dynamic> mat ;
    typedef Eigen::Matrix< std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> mat_complex ;
 
    class VectorFit
    {
    public:
        EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
        VectorFit();
        int Fitoption() const { return m_Fitoption; }
        void Fitoption(int val) { m_Fitoption = val; }
    
        bool RelaxedConstraint() const { return m_RelaxedConstraint; }
        void RelaxedConstraint(bool val) { m_RelaxedConstraint = val; }
        bool FlipUnstable() const { return m_FlipUnstable; }
        void FlipUnstable(bool val) { m_FlipUnstable = val; }
        bool PoleIdentification() const { return m_PoleIdentification; }
        void PoleIdentification(bool val) { m_PoleIdentification = val; }
        bool ResidueIdentification() const { return m_ResidueIdentification; }
        void ResidueIdentification(bool val) { m_ResidueIdentification = val; }
        bool ComplexStateSpace() const { return m_ComplexStateSpace; }
        void ComplexStateSpace(bool val) { m_ComplexStateSpace = val; }
    
        mat_complex MatrixA() const { return matrixA; }
        void MatrixA(mat_complex val) { matrixA = val; }
        mat_complex MatrixB() const { return matrixB; }
        void MatrixB(mat_complex val) { matrixB = val; }
        mat_complex MatrixC() const { return matrixC; }
        void MatrixC(mat_complex val) { matrixC = val; }
        std::complex<double> MatrixD() const { return matrixD; }
        void MatrixD(std::complex<double> val) { matrixD = val; }
        std::complex<double> MatrixE() const { return matrixE; }
        void MatrixE(std::complex<double> val) { matrixE = val; }
    
        mat MatrixAreal() const ;
        mat MatrixBreal() const ;
        mat MatrixCreal() const ;
        double MatrixDreal() const ;
        double MatrixEreal() const ;
    
        mat_complex Diff() const { return diff; }
        double Rmserr() const { return rmserr;}
    
        void Fit(const vec_complex &f, const vec_complex &s, vec_complex &poles, const vec & weight);
    
        vec_complex LinSolveLS( mat &A, mat &B);
        vec_complex LinSolve( mat &A, mat &B);
        vec_complex MatrixPolynomial( const mat_complex &A );
        vec_complex EigenValues( const mat_complex &ZER  );
    
        void QRFactor( mat &A_new, mat &Q, mat &R );    
        void TransferFunction( vec_complex& num, vec_complex &den);
    
    private:
        int     m_Fitoption;                // 1: Fitting with D=0,  E=0, 2: Fitting with D~=0, E=0, 3: Fitting with D~=0, E~=0 
        bool    m_RelaxedConstraint;        // Use relaxed nontriviality constraint 
        bool    m_FlipUnstable;             // Unstable poles are kept unchanged
        bool    m_PoleIdentification;       // The pole identification part is skipped, i.e (C,D,E) 
        bool    m_ResidueIdentification;    // The residue identification part is skipped, i.e. only the poles (A) are identified while C,D,E are returned as zero. 
        bool    m_ComplexStateSpace;        // The returned state-space model has real and complex  parameters. Output variable A is diagonal (and sparse). 
        
        mat_complex diff;
        double rmserr;
    
        mat_complex matrixA;
        mat_complex matrixB;
        mat_complex matrixC;
        std::complex<double> matrixD;
        std::complex<double> matrixE;
        static bool root_less_pred( std::complex<double> const& lh, std::complex<double> const &rh)
        {
            if( lh.real() != rh.real() )
                return lh.real() < rh.real();
            return lh.imag() < rh.imag();
        }
    
    };
}
 
 
#endif // CVectorfit_h__