ROOT logo
// @(#)root/mathmore:$Id: GSLNLSMinimizer.h 31604 2009-12-07 19:04:33Z moneta $
// Author: L. Moneta Wed Dec 20 17:16:32 2006

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 * This library is free software; you can redistribute it and/or      *
 * modify it under the terms of the GNU General Public License        *
 * as published by the Free Software Foundation; either version 2     *
 * of the License, or (at your option) any later version.             *
 *                                                                    *
 * This library is distributed in the hope that it will be useful,    *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
 * General Public License for more details.                           *
 *                                                                    *
 * You should have received a copy of the GNU General Public License  *
 * along with this library (see file COPYING); if not, write          *
 * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
 * 330, Boston, MA 02111-1307 USA, or contact the author.             *
 *                                                                    *
 **********************************************************************/

// Header file for class GSLNLSMinimizer

#ifndef ROOT_Math_GSLNLSMinimizer
#define ROOT_Math_GSLNLSMinimizer



#ifndef ROOT_Math_Minimizer
#include "Math/Minimizer.h"
#endif


#ifndef ROOT_Math_IFunctionfwd
#include "Math/IFunctionfwd.h"
#endif

#ifndef ROOT_Math_IParamFunctionfwd
#include "Math/IParamFunctionfwd.h"
#endif

#ifndef ROOT_Math_FitMethodFunction
#include "Math/FitMethodFunction.h"
#endif

#ifndef ROOT_Math_MinimizerVariable
#include "Math/MinimizerVariable.h"
#endif


#include <vector>
#include <map>
#include <string>

namespace ROOT { 

   namespace Math { 

      class GSLMultiFit; 


//________________________________________________________________________________
/** 
    LSResidualFunc class description. 
    Internal class used for accessing the residuals of the Least Square function
    and their derivates which are estimated numerically using GSL numerical derivation. 
    The class contains a pointer to the fit method function and an index specifying 
    the i-th residual and wraps it in a multi-dim gradient function interface
    ROOT::Math::IGradientFunctionMultiDim. 
    The class is used by ROOT::Math::GSLNLSMinimizer (GSL non linear least square fitter)

    @ingroup MultiMin
*/
class LSResidualFunc : public IMultiGradFunction { 
public: 

   //default ctor (required by CINT) 
   LSResidualFunc() : fIndex(0), fChi2(0)
   {}


   LSResidualFunc(const ROOT::Math::FitMethodFunction & func, unsigned int i) : 
      fIndex(i), 
      fChi2(&func), 
      fX2(std::vector<double>(func.NDim() ) )
   {}


   // copy ctor
   LSResidualFunc(const LSResidualFunc & rhs) :
      IMultiGenFunction(), 
      IMultiGradFunction() 
   { 
      operator=(rhs);
   } 

   // assignment
   LSResidualFunc & operator= (const LSResidualFunc & rhs) 
   { 
      fIndex = rhs.fIndex;
      fChi2 = rhs.fChi2; 
      fX2 = rhs.fX2;
      return *this;
   } 

   IMultiGenFunction * Clone() const { 
      return new LSResidualFunc(*fChi2,fIndex); 
   }

   unsigned int NDim() const { return fChi2->NDim(); }

   void Gradient( const double * x, double * g) const { 
      unsigned int n = NDim(); 
      std::copy(x,x+n,fX2.begin());
      const double kEps = 1.0E-4;
      double f0 = DoEval(x); 
      for (unsigned int i = 0; i < n; ++i) { 
         fX2[i] += kEps;
         g[i] =  ( DoEval(&fX2.front()) - f0 )/kEps;
         fX2[i] = x[i];
      }
   } 

private: 

   double DoEval (const double * x) const { 
      return fChi2->DataElement(x, fIndex);
   }
   
   double DoDerivative(const double * x, unsigned int icoord) const { 
      //return  ROOT::Math::Derivator::Eval(*this, x, icoord, 1E-8);
      std::copy(x,x+NDim(),fX2.begin());
      const double kEps = 1.0E-4;
      fX2[icoord] += kEps;
      return ( DoEval(&fX2.front()) - DoEval(x) )/kEps;
   }

   unsigned int fIndex; 
   const ROOT::Math::FitMethodFunction * fChi2; 
   mutable std::vector<double> fX2;  // cached vector
};


//_____________________________________________________________________________________________________
/** 
   GSLNLSMinimizer class for Non Linear Least Square fitting
   It Uses the Levemberg-Marquardt algorithm from 
   <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html">
   GSL Non Linear Least Square fitting</A>.

   @ingroup MultiMin
*/ 
class GSLNLSMinimizer : public  ROOT::Math::Minimizer {

public: 

   /** 
      Default constructor
   */ 
   GSLNLSMinimizer (int type = 0); 

   /** 
      Destructor (no operations)
   */ 
   ~GSLNLSMinimizer ();  

private:
   // usually copying is non trivial, so we make this unaccessible

   /** 
      Copy constructor
   */ 
   GSLNLSMinimizer(const GSLNLSMinimizer &) : ROOT::Math::Minimizer() {} 

   /** 
      Assignment operator
   */ 
   GSLNLSMinimizer & operator = (const GSLNLSMinimizer & rhs)  {
      if (this == &rhs) return *this;  // time saving self-test
      return *this;
   }

public: 

   /// set the function to minimize
   virtual void SetFunction(const ROOT::Math::IMultiGenFunction & func); 

   /// set gradient the function to minimize
   virtual void SetFunction(const ROOT::Math::IMultiGradFunction & func); 

   /// set free variable 
   virtual bool SetVariable(unsigned int ivar, const std::string & name, double val, double step); 

   /// set lower limited variable 
   virtual bool SetLowerLimitedVariable(unsigned int  ivar , const std::string & name , double val , double step , double lower );
   /// set upper limited variable
   virtual bool SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ); 
   /// set upper/lower limited variable 
   virtual bool SetLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower , double upper ); 
   /// set fixed variable
   virtual bool SetFixedVariable(unsigned int ivar , const std::string & name , double val );  

   /// set the value of an existing variable 
   virtual bool SetVariableValue(unsigned int ivar, double val );
   /// set the values of all existing variables (array must be dimensioned to the size of existing parameters)
   virtual bool SetVariableValues(const double * x);

   /// method to perform the minimization
   virtual  bool Minimize(); 

   /// return minimum function value
   virtual double MinValue() const { return fMinVal; } 

   /// return expected distance reached from the minimum
   virtual double Edm() const { return fEdm; } // not impl. }

   /// return  pointer to X values at the minimum 
   virtual const double *  X() const { return &fValues.front(); } 

   /// return pointer to gradient values at the minimum 
   virtual const double *  MinGradient() const; 

   /// number of function calls to reach the minimum 
   virtual unsigned int NCalls() const { return 0; } // not yet ipl.  

   /// this is <= Function().NDim() which is the total 
   /// number of variables (free+ constrained ones) 
   virtual unsigned int NDim() const { return fDim; }   

   /// number of free variables (real dimension of the problem) 
   /// this is <= Function().NDim() which is the total 
   virtual unsigned int NFree() const { return fNFree; }  

   /// minimizer provides error and error matrix
   virtual bool ProvidesError() const { return true; } 

   /// return errors at the minimum 
   virtual const double * Errors() const { return (fErrors.size() > 0) ? &fErrors.front() : 0; }
//  { 
//       static std::vector<double> err; 
//       err.resize(fDim);
//       return &err.front(); 
//    }

   /** return covariance matrices elements 
       if the variable is fixed the matrix is zero
       The ordering of the variables is the same as in errors
   */ 
   virtual double CovMatrix(unsigned int , unsigned int ) const;

   /// return covariance matrix status
   virtual int CovMatrixStatus() const;

protected: 

private: 
   

   unsigned int fDim;        // dimension of the function to be minimized 
   unsigned int fNFree;      // dimension of the internal function to be minimized 
   unsigned int fSize;        // number of fit points (residuals)
 

   ROOT::Math::GSLMultiFit * fGSLMultiFit;        // pointer to GSL multi fit solver 
   const ROOT::Math::FitMethodFunction * fObjFunc;      // pointer to Least square function
   
   double fMinVal;                                // minimum function value
   double fEdm;                                   // edm value
   double fLSTolerance;                           // Line Search Tolerance
   std::vector<double> fValues;           
   std::vector<double> fErrors;
   std::vector<double> fCovMatrix;              //  cov matrix (stored as cov[ i * dim + j] 
   std::vector<double> fSteps;
   std::vector<std::string> fNames;
   std::vector<LSResidualFunc> fResiduals;   //! transient Vector of the residual functions

   std::vector<ROOT::Math::EMinimVariableType> fVarTypes;  // vector specifyng the type of variables
   std::map< unsigned int, std::pair<double, double> > fBounds; // map specifying the bound using as key the parameter index


}; 

   } // end namespace Math

} // end namespace ROOT


#endif /* ROOT_Math_GSLNLSMinimizer */
 GSLNLSMinimizer.h:1
 GSLNLSMinimizer.h:2
 GSLNLSMinimizer.h:3
 GSLNLSMinimizer.h:4
 GSLNLSMinimizer.h:5
 GSLNLSMinimizer.h:6
 GSLNLSMinimizer.h:7
 GSLNLSMinimizer.h:8
 GSLNLSMinimizer.h:9
 GSLNLSMinimizer.h:10
 GSLNLSMinimizer.h:11
 GSLNLSMinimizer.h:12
 GSLNLSMinimizer.h:13
 GSLNLSMinimizer.h:14
 GSLNLSMinimizer.h:15
 GSLNLSMinimizer.h:16
 GSLNLSMinimizer.h:17
 GSLNLSMinimizer.h:18
 GSLNLSMinimizer.h:19
 GSLNLSMinimizer.h:20
 GSLNLSMinimizer.h:21
 GSLNLSMinimizer.h:22
 GSLNLSMinimizer.h:23
 GSLNLSMinimizer.h:24
 GSLNLSMinimizer.h:25
 GSLNLSMinimizer.h:26
 GSLNLSMinimizer.h:27
 GSLNLSMinimizer.h:28
 GSLNLSMinimizer.h:29
 GSLNLSMinimizer.h:30
 GSLNLSMinimizer.h:31
 GSLNLSMinimizer.h:32
 GSLNLSMinimizer.h:33
 GSLNLSMinimizer.h:34
 GSLNLSMinimizer.h:35
 GSLNLSMinimizer.h:36
 GSLNLSMinimizer.h:37
 GSLNLSMinimizer.h:38
 GSLNLSMinimizer.h:39
 GSLNLSMinimizer.h:40
 GSLNLSMinimizer.h:41
 GSLNLSMinimizer.h:42
 GSLNLSMinimizer.h:43
 GSLNLSMinimizer.h:44
 GSLNLSMinimizer.h:45
 GSLNLSMinimizer.h:46
 GSLNLSMinimizer.h:47
 GSLNLSMinimizer.h:48
 GSLNLSMinimizer.h:49
 GSLNLSMinimizer.h:50
 GSLNLSMinimizer.h:51
 GSLNLSMinimizer.h:52
 GSLNLSMinimizer.h:53
 GSLNLSMinimizer.h:54
 GSLNLSMinimizer.h:55
 GSLNLSMinimizer.h:56
 GSLNLSMinimizer.h:57
 GSLNLSMinimizer.h:58
 GSLNLSMinimizer.h:59
 GSLNLSMinimizer.h:60
 GSLNLSMinimizer.h:61
 GSLNLSMinimizer.h:62
 GSLNLSMinimizer.h:63
 GSLNLSMinimizer.h:64
 GSLNLSMinimizer.h:65
 GSLNLSMinimizer.h:66
 GSLNLSMinimizer.h:67
 GSLNLSMinimizer.h:68
 GSLNLSMinimizer.h:69
 GSLNLSMinimizer.h:70
 GSLNLSMinimizer.h:71
 GSLNLSMinimizer.h:72
 GSLNLSMinimizer.h:73
 GSLNLSMinimizer.h:74
 GSLNLSMinimizer.h:75
 GSLNLSMinimizer.h:76
 GSLNLSMinimizer.h:77
 GSLNLSMinimizer.h:78
 GSLNLSMinimizer.h:79
 GSLNLSMinimizer.h:80
 GSLNLSMinimizer.h:81
 GSLNLSMinimizer.h:82
 GSLNLSMinimizer.h:83
 GSLNLSMinimizer.h:84
 GSLNLSMinimizer.h:85
 GSLNLSMinimizer.h:86
 GSLNLSMinimizer.h:87
 GSLNLSMinimizer.h:88
 GSLNLSMinimizer.h:89
 GSLNLSMinimizer.h:90
 GSLNLSMinimizer.h:91
 GSLNLSMinimizer.h:92
 GSLNLSMinimizer.h:93
 GSLNLSMinimizer.h:94
 GSLNLSMinimizer.h:95
 GSLNLSMinimizer.h:96
 GSLNLSMinimizer.h:97
 GSLNLSMinimizer.h:98
 GSLNLSMinimizer.h:99
 GSLNLSMinimizer.h:100
 GSLNLSMinimizer.h:101
 GSLNLSMinimizer.h:102
 GSLNLSMinimizer.h:103
 GSLNLSMinimizer.h:104
 GSLNLSMinimizer.h:105
 GSLNLSMinimizer.h:106
 GSLNLSMinimizer.h:107
 GSLNLSMinimizer.h:108
 GSLNLSMinimizer.h:109
 GSLNLSMinimizer.h:110
 GSLNLSMinimizer.h:111
 GSLNLSMinimizer.h:112
 GSLNLSMinimizer.h:113
 GSLNLSMinimizer.h:114
 GSLNLSMinimizer.h:115
 GSLNLSMinimizer.h:116
 GSLNLSMinimizer.h:117
 GSLNLSMinimizer.h:118
 GSLNLSMinimizer.h:119
 GSLNLSMinimizer.h:120
 GSLNLSMinimizer.h:121
 GSLNLSMinimizer.h:122
 GSLNLSMinimizer.h:123
 GSLNLSMinimizer.h:124
 GSLNLSMinimizer.h:125
 GSLNLSMinimizer.h:126
 GSLNLSMinimizer.h:127
 GSLNLSMinimizer.h:128
 GSLNLSMinimizer.h:129
 GSLNLSMinimizer.h:130
 GSLNLSMinimizer.h:131
 GSLNLSMinimizer.h:132
 GSLNLSMinimizer.h:133
 GSLNLSMinimizer.h:134
 GSLNLSMinimizer.h:135
 GSLNLSMinimizer.h:136
 GSLNLSMinimizer.h:137
 GSLNLSMinimizer.h:138
 GSLNLSMinimizer.h:139
 GSLNLSMinimizer.h:140
 GSLNLSMinimizer.h:141
 GSLNLSMinimizer.h:142
 GSLNLSMinimizer.h:143
 GSLNLSMinimizer.h:144
 GSLNLSMinimizer.h:145
 GSLNLSMinimizer.h:146
 GSLNLSMinimizer.h:147
 GSLNLSMinimizer.h:148
 GSLNLSMinimizer.h:149
 GSLNLSMinimizer.h:150
 GSLNLSMinimizer.h:151
 GSLNLSMinimizer.h:152
 GSLNLSMinimizer.h:153
 GSLNLSMinimizer.h:154
 GSLNLSMinimizer.h:155
 GSLNLSMinimizer.h:156
 GSLNLSMinimizer.h:157
 GSLNLSMinimizer.h:158
 GSLNLSMinimizer.h:159
 GSLNLSMinimizer.h:160
 GSLNLSMinimizer.h:161
 GSLNLSMinimizer.h:162
 GSLNLSMinimizer.h:163
 GSLNLSMinimizer.h:164
 GSLNLSMinimizer.h:165
 GSLNLSMinimizer.h:166
 GSLNLSMinimizer.h:167
 GSLNLSMinimizer.h:168
 GSLNLSMinimizer.h:169
 GSLNLSMinimizer.h:170
 GSLNLSMinimizer.h:171
 GSLNLSMinimizer.h:172
 GSLNLSMinimizer.h:173
 GSLNLSMinimizer.h:174
 GSLNLSMinimizer.h:175
 GSLNLSMinimizer.h:176
 GSLNLSMinimizer.h:177
 GSLNLSMinimizer.h:178
 GSLNLSMinimizer.h:179
 GSLNLSMinimizer.h:180
 GSLNLSMinimizer.h:181
 GSLNLSMinimizer.h:182
 GSLNLSMinimizer.h:183
 GSLNLSMinimizer.h:184
 GSLNLSMinimizer.h:185
 GSLNLSMinimizer.h:186
 GSLNLSMinimizer.h:187
 GSLNLSMinimizer.h:188
 GSLNLSMinimizer.h:189
 GSLNLSMinimizer.h:190
 GSLNLSMinimizer.h:191
 GSLNLSMinimizer.h:192
 GSLNLSMinimizer.h:193
 GSLNLSMinimizer.h:194
 GSLNLSMinimizer.h:195
 GSLNLSMinimizer.h:196
 GSLNLSMinimizer.h:197
 GSLNLSMinimizer.h:198
 GSLNLSMinimizer.h:199
 GSLNLSMinimizer.h:200
 GSLNLSMinimizer.h:201
 GSLNLSMinimizer.h:202
 GSLNLSMinimizer.h:203
 GSLNLSMinimizer.h:204
 GSLNLSMinimizer.h:205
 GSLNLSMinimizer.h:206
 GSLNLSMinimizer.h:207
 GSLNLSMinimizer.h:208
 GSLNLSMinimizer.h:209
 GSLNLSMinimizer.h:210
 GSLNLSMinimizer.h:211
 GSLNLSMinimizer.h:212
 GSLNLSMinimizer.h:213
 GSLNLSMinimizer.h:214
 GSLNLSMinimizer.h:215
 GSLNLSMinimizer.h:216
 GSLNLSMinimizer.h:217
 GSLNLSMinimizer.h:218
 GSLNLSMinimizer.h:219
 GSLNLSMinimizer.h:220
 GSLNLSMinimizer.h:221
 GSLNLSMinimizer.h:222
 GSLNLSMinimizer.h:223
 GSLNLSMinimizer.h:224
 GSLNLSMinimizer.h:225
 GSLNLSMinimizer.h:226
 GSLNLSMinimizer.h:227
 GSLNLSMinimizer.h:228
 GSLNLSMinimizer.h:229
 GSLNLSMinimizer.h:230
 GSLNLSMinimizer.h:231
 GSLNLSMinimizer.h:232
 GSLNLSMinimizer.h:233
 GSLNLSMinimizer.h:234
 GSLNLSMinimizer.h:235
 GSLNLSMinimizer.h:236
 GSLNLSMinimizer.h:237
 GSLNLSMinimizer.h:238
 GSLNLSMinimizer.h:239
 GSLNLSMinimizer.h:240
 GSLNLSMinimizer.h:241
 GSLNLSMinimizer.h:242
 GSLNLSMinimizer.h:243
 GSLNLSMinimizer.h:244
 GSLNLSMinimizer.h:245
 GSLNLSMinimizer.h:246
 GSLNLSMinimizer.h:247
 GSLNLSMinimizer.h:248
 GSLNLSMinimizer.h:249
 GSLNLSMinimizer.h:250
 GSLNLSMinimizer.h:251
 GSLNLSMinimizer.h:252
 GSLNLSMinimizer.h:253
 GSLNLSMinimizer.h:254
 GSLNLSMinimizer.h:255
 GSLNLSMinimizer.h:256
 GSLNLSMinimizer.h:257
 GSLNLSMinimizer.h:258
 GSLNLSMinimizer.h:259
 GSLNLSMinimizer.h:260
 GSLNLSMinimizer.h:261
 GSLNLSMinimizer.h:262
 GSLNLSMinimizer.h:263
 GSLNLSMinimizer.h:264
 GSLNLSMinimizer.h:265
 GSLNLSMinimizer.h:266
 GSLNLSMinimizer.h:267
 GSLNLSMinimizer.h:268
 GSLNLSMinimizer.h:269
 GSLNLSMinimizer.h:270
 GSLNLSMinimizer.h:271
 GSLNLSMinimizer.h:272
 GSLNLSMinimizer.h:273
 GSLNLSMinimizer.h:274
 GSLNLSMinimizer.h:275
 GSLNLSMinimizer.h:276
 GSLNLSMinimizer.h:277
 GSLNLSMinimizer.h:278
 GSLNLSMinimizer.h:279
 GSLNLSMinimizer.h:280
 GSLNLSMinimizer.h:281
 GSLNLSMinimizer.h:282
 GSLNLSMinimizer.h:283
 GSLNLSMinimizer.h:284
 GSLNLSMinimizer.h:285
 GSLNLSMinimizer.h:286
 GSLNLSMinimizer.h:287
 GSLNLSMinimizer.h:288
 GSLNLSMinimizer.h:289
 GSLNLSMinimizer.h:290
 GSLNLSMinimizer.h:291