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>
#include <ICLUtils/SmartPtr.h>
#include <ICLUtils/Function.h>

namespace icl{
  namespace math{


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

      typedef utils::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;

      utils::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 &){
          Si = m_C ? m_S.pinv(true)* (*m_C) : m_S.pinv(true);
        }


        // might cause an exception
        Si.eigen(m_Evecs, m_Evals);
        std::copy(m_Evecs.col_begin(0), m_Evecs.col_end(0), m_model.begin());

        return m_model;
      }
    };
  } // namespace math
}