LeastSquareModelFitting.h

Go to the documentation of this file.

/********************************************************************
**                Image Component Library (ICL)                    **
**                                                                 **
** Copyright (C) 2006-2013 CITEC, University of Bielefeld          **
**                         Neuroinformatics Group                  **
** Website: www.iclcv.org and                                      **
**          http://opensource.cit-ec.de/projects/icl               **
**                                                                 **
** File   : ICLMath/src/ICLMath/LeastSquareModelFitting.h          **
** Module : ICLMath                                                **
** Authors: Christof Elbrechter                                    **
**                                                                 **
**                                                                 **
** GNU LESSER GENERAL PUBLIC LICENSE                               **
** This file may be used under the terms of the GNU Lesser General **
** Public License version 3.0 as published by the                  **
**                                                                 **
** Free Software Foundation and appearing in the file LICENSE.LGPL **
** included in the packaging of this file.  Please review the      **
** following information to ensure the license requirements will   **
** be met: http://www.gnu.org/licenses/lgpl-3.0.txt                **
**                                                                 **
** The development of this software was supported by the           **
** Excellence Cluster EXC 277 Cognitive Interaction Technology.    **
** The Excellence Cluster EXC 277 is a grant of the Deutsche       **
** Forschungsgemeinschaft (DFG) in the context of the German       **
** Excellence Initiative.                                          **
**                                                                 **
********************************************************************/

#pragma once

#include <ICLUtils/CompatMacros.h>
#include <ICLMath/DynMatrix.h>

namespace icl{
  namespace math{
    

    template<class T, class DataPoint>
    class LeastSquareModelFitting{
      public:

      typedef Function<void,const DataPoint&,T*> DesignMatrixGen;
      
      typedef std::vector<T> Model;
      
      private:
      int m_modelDim;
  
      DesignMatrixGen m_gen;
      
      DynMatrix<T> m_D, m_S, m_Evecs, m_Evals, m_svdU, m_svdS, m_svdVt;
      
      SmartPtr<DynMatrix<T> >m_C;
      
      Model m_model;
      
      public:
      
      LeastSquareModelFitting(){}
      
      LeastSquareModelFitting(int modelDim, DesignMatrixGen gen, 
                              DynMatrix<T> *constraintMatrix=0):
      m_modelDim(modelDim),m_gen(gen),m_S(modelDim,modelDim), 
      m_C(constraintMatrix),m_model(modelDim){
        
      }
      

      icl64f getError(const Model &model,const DataPoint &p){
        std::vector<T> d(m_modelDim);
        m_gen(p,d.data());
        icl64f e = 0;
        for(int i=0;i<m_modelDim;++i) e += d[i] * model[i];
        return fabs(e);
      }
      

      Model fit(const std::vector<DataPoint> &points){
        const int M = m_modelDim;
        const int N = (int)points.size();
        
        m_D.setBounds(M,N);
  
        for(int i=0;i<N;i++){
          m_gen(points[i],m_D.row_begin(i));
        }
        
        m_D.transp().mult(m_D,m_S);
        
        
        
        DynMatrix<T> Si;
        try{ 
          Si = m_C ? m_S.inv()* (*m_C) : m_S.inv();
        }catch(SingularMatrixException &ex){
          Si = m_C ? m_S.pinv(true)* (*m_C) : m_S.pinv(true); 
        }
        
        try{
          Si.eigen(m_Evecs, m_Evals);
          std::copy(m_Evecs.col_begin(0), m_Evecs.col_end(0), m_model.begin());
        }catch(ICLException &e){
          //Si.svd(m_svdU, m_svdS, m_svdVt);
          //SHOW(m_svdS);
          //
          //std::copy(m_svdU.col_begin(0), m_svdU.col_end(0), m_model.begin());
          //
          //for(int i=0;i<N;++i){
          //  SHOW(getError(m_model,points[i]));
          //}

          // this can happen if there is no noise in the input values
          // rather than solving Ax = lambda x, we have a 0 eigenvalue lambda
          // therefore we have to solve Ax = 0 ??
          DEBUG_LOG("LeastSquareModelFitting:fit: error in eigenvalue decomposition:\n" << e.what());
          std::fill(m_model.begin(),m_model.end(),Range<T>::limits().maxVal);
        }
        
        return m_model;
      }
    };
  } // namespace math
}