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