FumiliGradientCalculator.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id$
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
 *                                                                    *
 **********************************************************************/
 
#include "Minuit2/FumiliGradientCalculator.h"
#include "Minuit2/FumiliFCNBase.h"
#include "Minuit2/MnUserTransformation.h"
#include "Minuit2/FunctionGradient.h"
#include "Minuit2/MinimumParameters.h"
#include "Minuit2/FumiliChi2FCN.h"
#include "Minuit2/FumiliMaximumLikelihoodFCN.h"
 
//to compare with N2P calculator
//#define DEBUG 1
#ifdef DEBUG
#include "Minuit2/MnPrint.h"
#include "Minuit2/Numerical2PGradientCalculator.h"
#include "Minuit2/MnStrategy.h"
#include "Minuit2/MnUserFcn.h"
#endif
 
namespace ROOT {
 
   namespace Minuit2 {
 
 
FunctionGradient FumiliGradientCalculator::operator()(const MinimumParameters& par) const {
 
   // Calculate gradient for Fumili using the gradient and Hessian provided by the FCN Fumili function
   // applying the external to int trasformation.
 
   int nvar = par.Vec().size();
   std::vector<double> extParam = fTransformation(  par.Vec() );
   //   std::vector<double> deriv;
   //   deriv.reserve( nvar );
   //   for (int i = 0; i < nvar; ++i) {
   //     unsigned int ext = fTransformation.ExtOfInt(i);
   //     if ( fTransformation.Parameter(ext).HasLimits())
   //       deriv.push_back( fTransformation.DInt2Ext( i, par.Vec()(i) ) );
   //     else
   //       deriv.push_back(1.0);
   //   }
 
   // eval Gradient
   FumiliFCNBase & fcn = const_cast<FumiliFCNBase &>(fFcn);
 
   fcn.EvaluateAll(extParam);
 
 
   MnAlgebraicVector v(nvar);
   MnAlgebraicSymMatrix h(nvar);
 
 
   const std::vector<double> & fcn_gradient = fFcn.Gradient();
   assert( fcn_gradient.size() == extParam.size() );
 
 
   //   for (int i = 0; i < nvar; ++i) {
   //     unsigned int iext = fTransformation.ExtOfInt(i);
   //     double ideriv = 1.0;
   //     if ( fTransformation.Parameter(iext).HasLimits())
   //       ideriv =  fTransformation.DInt2Ext( i, par.Vec()(i) ) ;
 
 
   //     //     v(i) = fcn_gradient[iext]*deriv;
   //     v(i) = ideriv*fcn_gradient[iext];
 
   //     for (int j = i; j < nvar; ++j) {
   //       unsigned int jext = fTransformation.ExtOfInt(j);
   //       double jderiv = 1.0;
   //       if ( fTransformation.Parameter(jext).HasLimits())
   // jderiv =  fTransformation.DInt2Ext( j, par.Vec()(j) ) ;
 
   // //       h(i,j) = deriv[i]*deriv[j]*fFcn.Hessian(iext,jext);
   //       h(i,j) = ideriv*jderiv*fFcn.Hessian(iext,jext);
   //     }
   //   }
 
 
   // cache deriv and Index values .
   // in large Parameter limit then need to re-optimize and see if better not caching
 
   std::vector<double> deriv(nvar);
   std::vector<unsigned int> extIndex(nvar);
   for (int i = 0; i < nvar; ++i) {
      extIndex[i] = fTransformation.ExtOfInt(i);
      deriv[i] = 1;
      if ( fTransformation.Parameter(extIndex[i]).HasLimits())
         deriv[i] =  fTransformation.DInt2Ext( i, par.Vec()(i) ) ;
 
      v(i) = fcn_gradient[extIndex[i]]*deriv[i];
 
      for (int j = 0; j <= i; ++j) {
         h(i,j) = deriv[i]*deriv[j]*fFcn.Hessian(extIndex[i],extIndex[j]);
      }
   }
 
#ifdef DEBUG
   // compare with other Gradient
   //   // calculate Gradient using Minuit method
 
   Numerical2PGradientCalculator gc(MnUserFcn(fFcn,fTransformation), fTransformation, MnStrategy(1));
   FunctionGradient g2 = gc(par);
 
   std::cout << "Fumili Gradient " << v << std::endl;
   std::cout << "Minuit Gradient " << g2.Vec() << std::endl;
#endif
 
   // store calculated Hessian
   fHessian = h;
   return FunctionGradient(v);
}
 
FunctionGradient FumiliGradientCalculator::operator()(const MinimumParameters& par,
                                                      const FunctionGradient&) const
 
{
   // Needed for interface of base class.
   return this->operator()(par);
 
}
 
   }  // namespace Minuit2
 
}  // namespace ROOT