LevenbergMarquardtFitter.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/LevenbergMarquardtFitter.h         **
** Module : ICLMath                                                **
** Authors: Christof Elbrechter, Sergius Gaulik                    **
**                                                                 **
**                                                                 **
** 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/DynVector.h>
#include <ICLUtils/Function.h>

namespace icl{
  namespace math{



    template<class Scalar>
    class ICLMath_API LevenbergMarquardtFitter{
      public:

      typedef DynColVector<Scalar> Vector; 
      typedef DynColVector<Scalar> Params; 
      typedef DynMatrix<Scalar> Matrix;    

      struct Result{
        int iteration;                
        Scalar initialError;          
        Scalar error;                 
        std::vector<Scalar> lambdas;  
        Params params;                

        friend ICLMath_API inline std::ostream &operator<<(std::ostream &str, const Result &d){
          return str << "iteration: " << d.iteration << " initial error: " << d.initialError << " error:" << d.error
                     << "  lambda[0]:" << d.lambdas[0] << "  params:" << d.params.transp();
        }
      };

      struct Data{
        Matrix x;  
        Matrix y;  
      };

      typedef icl::utils::Function<Vector,const Params&,const Vector&> Function;
      typedef icl::utils::Function<Matrix,const Params&,const Matrix&> FunctionMat;


      typedef icl::utils::Function<void,const Params&, const Vector&, Vector&> Jacobian;
      typedef icl::utils::Function<void,const Params&, const Matrix&, Matrix&> JacobianMat;

      typedef icl::utils::Function<void,const Result&> DebugCallback;



      protected:
      Function f;        
      FunctionMat fMat;  

      bool useMultiThreading;   
      bool useMat;              
      Scalar tau;               
      int maxIterations;        
      Scalar minError;          
      Scalar lambdaMultiplier;  
      Scalar eps1;              
      Scalar eps2    ;          
      std::string linSolver;    

      Vector dst;        
      Matrix J;          
      Matrix H;          

      std::vector<Jacobian> js;
      std::vector<JacobianMat> jsMat;

      DebugCallback dbg;  
      Params params_new;  
      Matrix y_est;       
      //Matrix y_est_tmp;   //!< current estimated outputs (with temporary update)
      Vector dy;          


      Scalar error(const Matrix &ys, const Matrix &y_est) const;

      public:

      LevenbergMarquardtFitter();


      LevenbergMarquardtFitter(Function f, int outputDim,
          const std::vector<Jacobian> &js=std::vector<Jacobian>(),
          Scalar tau=1.e-3, int maxIterations=200,
          Scalar minError = 1.e-6, Scalar lambdaMultiplier=10,
          Scalar eps1 = 1.49012e-08, Scalar eps2 = 1.49012e-08,
          const std::string &linSolver="svd");
      LevenbergMarquardtFitter(FunctionMat f, int outputDim,
          const std::vector<JacobianMat> &js=std::vector<JacobianMat>(),
          Scalar tau=1.e-3, int maxIterations=200,
          Scalar minError = 1.e-6, Scalar lambdaMultiplier=10,
          Scalar eps1 = 1.49012e-08, Scalar eps2 = 1.49012e-08,
          const std::string &linSolver="svd");


      void setUseMultiThreading(bool enable);


      void init(Function f, int outputDim,
                const std::vector<Jacobian> &js=std::vector<Jacobian>(),
                Scalar tau=1.e-8, int maxIterations=1000,
                Scalar minError = 1.e-6, Scalar lambdaMultiplier=10,
                Scalar eps1 = 1.49012e-08, Scalar eps2 = 1.49012e-08,
                const std::string &linSolver="svd");
      void init(FunctionMat f, int outputDim,
                const std::vector<JacobianMat> &js=std::vector<JacobianMat>(),
                Scalar tau=1.e-8, int maxIterations=1000,
                Scalar minError = 1.e-6, Scalar lambdaMultiplier=10,
                Scalar eps1 = 1.49012e-08, Scalar eps2 = 1.49012e-08,
                const std::string &linSolver="svd");


      Result fit(const Matrix &xs, const Matrix &ys, Params initParams);




      static Jacobian create_numerical_jacobian(int o, Function f, float delta=1.E-5);
      static JacobianMat create_numerical_jacobian(int o, FunctionMat f, float delta=1.E-5);

      static std::vector<Jacobian> create_numerical_jacobians(int n, Function f, float delta=1.e-5);
      static std::vector<JacobianMat> create_numerical_jacobians(int n, FunctionMat f, float delta=1.e-5);


      static Data create_data(const Params &p, Function f, int xDim, int yDim,
                              int num=1000, Scalar minX=-5, Scalar maxX=5);

      static void default_debug_callback(const Result &r);

      void setDebugCallback(DebugCallback dbg=default_debug_callback);

    private:
      Result fitVec(const Matrix &xs, const Matrix &ys, Params initParams);
      Result fitMat(const Matrix &xs, const Matrix &ys, Params initParams);
    };
  } // namespace math
} // namespace icl