Galaxis Codebase: PolyFit< T > Class Template Reference

#include <polyfit.h>

Collaboration diagram for PolyFit< T >:

Classes
struct	Options

Public Types
enum	SolverType { EIGEN_COL_PIV_HOUSEHOLDER, EIGEN_FULL_PIV_HOUSEHOLDER, EIGEN_JACOBI_SVD }

Public Member Functions
int	getDegree () const

SolverType	getSolverType () const

const VecType	getCoefficients () const

	PolyFit (std::vector< T > &x, std::vector< T > &y, Options &options)

std::vector< T >	evalPoly (const std::vector< T > &at)

void	setSamplePoints (std::vector< T > &v)

void	setSignalPoints (std::vector< T > &v)

T	solveLS (void)

T	getMSE (void)

void	randsample (int N, int s, std::set< int > &ind)

void	solveRLS (void)

Private Types
typedef boost::uniform_int< int >	DistributionType

typedef Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic >	MatType

typedef Eigen::Matrix< T, Eigen::Dynamic, 1 >	VecType

Private Member Functions
void	dataToEigenSubset (std::set< int > &indexSet)

int	getInliers (std::vector< int > &inliers)

void	dataToEigen (void)

void	buildV (void)

Private Attributes
boost::mt19937	gen

VecType	pCoef

std::vector< T > &	x

std::vector< T > &	y

VecType	xM

VecType	yM

MatType	V

Options &	options

Detailed Description

template<typename T>
class PolyFit< T >

Definition at line 38 of file polyfit.h.

Member Typedef Documentation

◆ DistributionType

template<typename T >

typedef boost::uniform_int<int> PolyFit< T >::DistributionType

private

Definition at line 69 of file polyfit.h.

◆ MatType

template<typename T >

typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> PolyFit< T >::MatType

private

Definition at line 73 of file polyfit.h.

◆ VecType

template<typename T >

typedef Eigen::Matrix<T, Eigen::Dynamic, 1> PolyFit< T >::VecType

private

Definition at line 76 of file polyfit.h.

Member Enumeration Documentation

◆ SolverType

template<typename T >

enum PolyFit::SolverType

Enumerator
EIGEN_COL_PIV_HOUSEHOLDER
EIGEN_FULL_PIV_HOUSEHOLDER
EIGEN_JACOBI_SVD

Definition at line 40 of file polyfit.h.

                   {
     EIGEN_COL_PIV_HOUSEHOLDER,
     EIGEN_FULL_PIV_HOUSEHOLDER,
     EIGEN_JACOBI_SVD
   };

Constructor & Destructor Documentation

◆ PolyFit()

template<typename T >

PolyFit< T >::PolyFit	(	std::vector< T > &	x,
		std::vector< T > &	y,
		Options &	options
	)

inline

Definition at line 109 of file polyfit.h.

                                                               :
       x(x), y(y), options(options) {
     if(!x.size() || !y.size() || (x.size() != y.size())) {
       throw std::exception();
     }
  
     if(options.polyDeg < 1) {
       throw std::exception();
     }
     dataToEigen();
   }

References PolyFit< T >::dataToEigen(), PolyFit< T >::options, PolyFit< T >::Options::polyDeg, PolyFit< T >::x, and PolyFit< T >::y.

Here is the call graph for this function:

Member Function Documentation

◆ buildV()

template<typename T >

void PolyFit< T >::buildV ( void )

inlineprivate

Definition at line 321 of file polyfit.h.

                     {
     V.resize(xM.rows(), options.polyDeg + 1);
  
     for(int i = 0; i < xM.rows(); ++i)
       V(i, options.polyDeg) = 1;
  
     for(int d = 0; d < options.polyDeg; ++d)
       for(int i = 0; i < xM.rows(); ++i)
         V(i, d) = pow(xM(i, 0), options.polyDeg - d);
   }

References PolyFit< T >::options, PolyFit< T >::Options::polyDeg, PolyFit< T >::V, and PolyFit< T >::xM.

Referenced by PolyFit< T >::getInliers(), PolyFit< T >::solveLS(), and PolyFit< T >::solveRLS().

Here is the caller graph for this function:

◆ dataToEigen()

template<typename T >

void PolyFit< T >::dataToEigen ( void )

inlineprivate

Definition at line 309 of file polyfit.h.

                          {
     xM.resize(static_cast<int>(x.size()), 1);
     yM.resize(static_cast<int>(y.size()), 1);
  
     for(int i = 0; i < static_cast<int>(x.size()); ++i)
       xM(i, 0) = x[i];
  
     for(int i = 0; i < static_cast<int>(y.size()); ++i)
       yM(i, 0) = y[i];
   }

References PolyFit< T >::x, PolyFit< T >::xM, PolyFit< T >::y, and PolyFit< T >::yM.

Referenced by PolyFit< T >::getInliers(), PolyFit< T >::PolyFit(), PolyFit< T >::setSamplePoints(), and PolyFit< T >::setSignalPoints().

Here is the caller graph for this function:

◆ dataToEigenSubset()

template<typename T >

void PolyFit< T >::dataToEigenSubset ( std::set< int > & indexSet )

inlineprivate

Definition at line 278 of file polyfit.h.

                                                 {
     int nVal = static_cast<int>(indexSet.size());
     xM.resize(nVal, 1);
     yM.resize(nVal, 1);
  
     int cnt                    = 0;
     std::set<int>::iterator it = indexSet.begin();
     while(it != indexSet.end()) {
       xM(cnt, 0) = x[*it];
       yM(cnt, 0) = y[*it];
       ++cnt;
       ++it;
     }
   }

References PolyFit< T >::x, PolyFit< T >::xM, PolyFit< T >::y, and PolyFit< T >::yM.

Referenced by PolyFit< T >::solveRLS().

Here is the caller graph for this function:

◆ evalPoly()

template<typename T >

std::vector<T> PolyFit< T >::evalPoly ( const std::vector< T > & at )

inline

Definition at line 122 of file polyfit.h.

                                                 {
     VecType x;
     VecType y;
  
     x.resize(at.size());
     y.resize(at.size());
  
     for(int i = 0; i < at.size(); ++i)
       x(i, 0) = at[i];
  
     // Horner method
     if(pCoef.rows() > 0)
       y.setConstant(pCoef(0, 0));
  
     int nc = pCoef.rows();
     for(int i = 1; i < nc; ++i)
       y = (x.cwiseProduct(y)).array() + pCoef(i, 0);
  
     // Copy data over to output vector
     std::vector<T> out(y.rows());
     for(int i = 0; i < y.rows(); ++i)
       out[i] = y(i, 0);
  
     return out;
   }

References PolyFit< T >::pCoef, PolyFit< T >::x, and PolyFit< T >::y.

◆ getCoefficients()

template<typename T >

const VecType PolyFit< T >::getCoefficients ( ) const

inline

Definition at line 104 of file polyfit.h.

                                         {
     return pCoef;
   }

References PolyFit< T >::pCoef.

Referenced by operator<<().

Here is the caller graph for this function:

◆ getDegree()

template<typename T >

int PolyFit< T >::getDegree ( ) const

inline

Definition at line 96 of file polyfit.h.

                         {
     return options.polyDeg;
   }

References PolyFit< T >::options, and PolyFit< T >::Options::polyDeg.

Referenced by operator<<().

Here is the caller graph for this function:

◆ getInliers()

template<typename T >

int PolyFit< T >::getInliers ( std::vector< int > & inliers )

inlineprivate

Definition at line 295 of file polyfit.h.

                                           {
     // Full data to Eigen
     dataToEigen();
     buildV();
  
     Eigen::VectorXd diff = ((V * pCoef).array() - yM.col(0).array());
     for(int i = 0; i < diff.rows(); ++i) {
       if(pow(diff(i), 2) < options.tolerance)
         inliers.push_back(i);
     }
     return static_cast<int>(inliers.size());
   }

References PolyFit< T >::buildV(), PolyFit< T >::dataToEigen(), PolyFit< T >::options, PolyFit< T >::pCoef, PolyFit< T >::Options::tolerance, PolyFit< T >::V, and PolyFit< T >::yM.

Referenced by PolyFit< T >::solveRLS().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getMSE()

template<typename T >

T PolyFit< T >::getMSE ( void )

inline

Definition at line 199 of file polyfit.h.

                  {
     if(V.cols() != pCoef.rows() && yM.rows() != V.rows()) {
       throw std::logic_error("Assert Polyfit l220");
     }
  
     T mse                = 0.0;
     Eigen::VectorXd diff = ((V * pCoef).array() - yM.col(0).array());
     for(int i = 0; i < diff.rows(); ++i)
       mse += pow(diff(i), 2);
     return mse;
   }

References PolyFit< T >::pCoef, PolyFit< T >::V, and PolyFit< T >::yM.

Referenced by PolyFit< T >::solveLS().

Here is the caller graph for this function:

◆ getSolverType()

template<typename T >

SolverType PolyFit< T >::getSolverType ( ) const

inline

Definition at line 100 of file polyfit.h.

                                    {
     return options.solver;
   }

References PolyFit< T >::options, and PolyFit< T >::Options::solver.

Referenced by operator<<().

Here is the caller graph for this function:

◆ randsample()

template<typename T >

void PolyFit< T >::randsample	(	int	N,
		int	s,
		std::set< int > &	ind
	)

inline

Definition at line 212 of file polyfit.h.

                                                   {
     if(N < 0 && s < 0 && (N <= s)) {
       throw std::logic_error("Assert Polyfit l230");
     }
  
     DistributionType dist(1, N);
     while(static_cast<int>(ind.size()) < s) {
       ind.insert(dist(gen));
     }
   }

References PolyFit< T >::gen.

Referenced by PolyFit< T >::solveRLS().

Here is the caller graph for this function:

◆ setSamplePoints()

template<typename T >

void PolyFit< T >::setSamplePoints ( std::vector< T > & v )

inline

Definition at line 149 of file polyfit.h.

                                         {
     if(!v.size()) {
       throw std::exception();
     }
  
     if(x.size()) {
     }
     x = v;
     dataToEigen();
   }

References PolyFit< T >::dataToEigen(), and PolyFit< T >::x.

Here is the call graph for this function:

◆ setSignalPoints()

template<typename T >

void PolyFit< T >::setSignalPoints ( std::vector< T > & v )

inline

Definition at line 161 of file polyfit.h.

                                         {
     if(!v.size()) {
       throw std::exception();
     }
  
     if(y.size()) {
     }
     y = v;
     dataToEigen();
   }

References PolyFit< T >::dataToEigen(), and PolyFit< T >::y.

Here is the call graph for this function:

◆ solveLS()

template<typename T >

T PolyFit< T >::solveLS ( void )

inline

Definition at line 173 of file polyfit.h.

                   {
  
     // Size check
     if(xM.rows() != yM.rows()) {
       throw std::exception();
     }
  
     buildV();
  
     switch(options.solver) {
       case EIGEN_COL_PIV_HOUSEHOLDER:
         pCoef = V.colPivHouseholderQr().solve(yM.col(0));
         break;
       case EIGEN_FULL_PIV_HOUSEHOLDER:
         pCoef = V.fullPivHouseholderQr().solve(yM.col(0));
         break;
       case EIGEN_JACOBI_SVD:
         pCoef = V.jacobiSvd(Eigen::ComputeThinU |
                             Eigen::ComputeThinV)
                     .solve(yM.col(0));
         break;
     }
     return getMSE();
   }

References PolyFit< T >::buildV(), PolyFit< T >::EIGEN_COL_PIV_HOUSEHOLDER, PolyFit< T >::EIGEN_FULL_PIV_HOUSEHOLDER, PolyFit< T >::EIGEN_JACOBI_SVD, PolyFit< T >::getMSE(), PolyFit< T >::options, PolyFit< T >::pCoef, PolyFit< T >::Options::solver, PolyFit< T >::V, PolyFit< T >::xM, and PolyFit< T >::yM.

Referenced by PolyFit< T >::solveRLS().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ solveRLS()

template<typename T >

void PolyFit< T >::solveRLS ( void )

inline

Definition at line 224 of file polyfit.h.

                       {
     std::vector<int> bestSubset;
     Eigen::Matrix<T, Eigen::Dynamic, 1> bestCoef;
  
     // Probability that we can at least choose
     // one sample with options.maxMPts free of
     // outliers (sugg. by Peter Kovesi's code)
     double p = 0.99;
  
     int trialCount = 0;
     int N          = 1;
     while(N > trialCount) {
  
       std::set<int> indexSet;
       randsample(static_cast<int>(y.size() - 1),
                  options.maxMPts, indexSet);
  
       if(static_cast<int>(indexSet.size()) != options.maxMPts) {
         throw std::logic_error("Assert Polyfit l265");
       }
  
       // Build (sample,signal) for current subset
       dataToEigenSubset(indexSet);
       buildV();
  
       // Solve the LS-fitting problem
       double mse = solveLS();
  
       // Compute inlier indices
       std::vector<int> inliers;
       int nInliers = getInliers(inliers);
  
       if(inliers.size() > bestSubset.size()) {
         // Store settings
         bestSubset = inliers;
         bestCoef   = pCoef;
  
         // Refine estimate of N (based on Peter Kovesi's MATLAB impl. of RANSAC)
         double fracInliers = static_cast<double>(inliers.size()) / double(yM.rows());
         double pNoOutliers = 1 - pow(fracInliers, options.maxMPts);
         pNoOutliers        = std::max(1e-9, pNoOutliers);
         pNoOutliers        = std::min(1 - 1e-9, pNoOutliers);
         N                  = log(1 - p) / log(pNoOutliers);
       }
  
       ++trialCount;
       if(trialCount > options.maxTrial) {
         break;
       }
     }
     pCoef = bestCoef;
   }