MnSeedGenerator.cxx

Go to the documentation of this file.
00001 // @(#)root/minuit2:$Id: MnSeedGenerator.cxx 23522 2008-04-24 15:09:19Z moneta $
00002 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005  
00003 
00004 /**********************************************************************
00005  *                                                                    *
00006  * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
00007  *                                                                    *
00008  **********************************************************************/
00009 
00010 #include "Minuit2/MnSeedGenerator.h"
00011 #include "Minuit2/MinimumSeed.h"
00012 #include "Minuit2/MnFcn.h"
00013 #include "Minuit2/GradientCalculator.h"
00014 #include "Minuit2/InitialGradientCalculator.h"
00015 #include "Minuit2/MnUserTransformation.h"
00016 #include "Minuit2/MinimumParameters.h"
00017 #include "Minuit2/FunctionGradient.h"
00018 #include "Minuit2/MinimumError.h"
00019 #include "Minuit2/MnMatrix.h"
00020 #include "Minuit2/MnMachinePrecision.h"
00021 #include "Minuit2/MinuitParameter.h"
00022 #include "Minuit2/MnLineSearch.h"
00023 #include "Minuit2/MnParabolaPoint.h"
00024 #include "Minuit2/MinimumState.h"
00025 #include "Minuit2/MnUserParameterState.h"
00026 #include "Minuit2/MnStrategy.h"
00027 #include "Minuit2/MnHesse.h"
00028 #include "Minuit2/VariableMetricEDMEstimator.h"
00029 #include "Minuit2/NegativeG2LineSearch.h"
00030 #include "Minuit2/AnalyticalGradientCalculator.h"
00031 #include "Minuit2/Numerical2PGradientCalculator.h"
00032 #include "Minuit2/HessianGradientCalculator.h"
00033 
00034 //#define DEBUG
00035 
00036 #if defined(DEBUG) || defined(WARNINGMSG)
00037 #include "Minuit2/MnPrint.h"
00038 #endif
00039 
00040 
00041 
00042 #include <math.h>
00043 
00044 
00045 namespace ROOT {
00046 
00047    namespace Minuit2 {
00048 
00049 
00050 MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculator& gc, const MnUserParameterState& st, const MnStrategy& stra) const {
00051 
00052 
00053    // find seed (initial minimization point) using the calculated gradient
00054    unsigned int n = st.VariableParameters();
00055    const MnMachinePrecision& prec = st.Precision();
00056 
00057 #ifdef DEBUG
00058    std::cout << "MnSeedGenerator: operator() - var par = " << n << " mnfcn pointer " << &fcn << std::endl;
00059 #endif
00060    
00061    // initial starting values
00062    MnAlgebraicVector x(n);
00063    for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
00064    double fcnmin = fcn(x);
00065    MinimumParameters pa(x, fcnmin);
00066    FunctionGradient dgrad = gc(pa);
00067    MnAlgebraicSymMatrix mat(n);
00068    double dcovar = 1.;
00069    if(st.HasCovariance()) {
00070       for(unsigned int i = 0; i < n; i++)       
00071          for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
00072       dcovar = 0.;
00073    } else {
00074       for(unsigned int i = 0; i < n; i++)       
00075          mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
00076    }
00077    MinimumError err(mat, dcovar);
00078    double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
00079    MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
00080    
00081    NegativeG2LineSearch ng2ls;
00082    if(ng2ls.HasNegativeG2(dgrad, prec)) {
00083 #ifdef DEBUG
00084       std::cout << "MnSeedGenerator: Negative G2 Found: " << std::endl;
00085       std::cout << x << std::endl; 
00086       std::cout << dgrad.Grad() << std::endl; 
00087       std::cout << dgrad.G2() << std::endl; 
00088 #endif
00089       state = ng2ls(fcn, state, gc, prec);
00090    }
00091 
00092    
00093    if(stra.Strategy() == 2 && !st.HasCovariance()) {
00094       //calculate full 2nd derivative
00095 #ifdef DEBUG
00096       std::cout << "MnSeedGenerator: calling MnHesse  " << std::endl;
00097 #endif
00098       MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
00099       return MinimumSeed(tmp, st.Trafo());
00100    }
00101    
00102    return MinimumSeed(state, st.Trafo());
00103 }
00104 
00105 
00106 MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradientCalculator& gc, const MnUserParameterState& st, const MnStrategy& stra) const {
00107    // find seed (initial point for minimization) using analytical gradient
00108    unsigned int n = st.VariableParameters();
00109    const MnMachinePrecision& prec = st.Precision();
00110    
00111    // initial starting values
00112    MnAlgebraicVector x(n);
00113    for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
00114    double fcnmin = fcn(x);
00115    MinimumParameters pa(x, fcnmin);
00116    
00117    InitialGradientCalculator igc(fcn, st.Trafo(), stra);
00118    FunctionGradient tmp = igc(pa);
00119    FunctionGradient grd = gc(pa);
00120    FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep());
00121    
00122    if(gc.CheckGradient()) {
00123       bool good = true;
00124       HessianGradientCalculator hgc(fcn, st.Trafo(), MnStrategy(2));
00125       std::pair<FunctionGradient, MnAlgebraicVector> hgrd = hgc.DeltaGradient(pa, dgrad);
00126       for(unsigned int i = 0; i < n; i++) {
00127          if(fabs(hgrd.first.Grad()(i) - grd.Grad()(i)) > hgrd.second(i)) {
00128 #ifdef WARNINGMSG
00129             MN_INFO_MSG("MnSeedGenerator:gradient discrepancy of external Parameter too large");
00130             int externalParameterIndex = st.Trafo().ExtOfInt(i); 
00131             MN_INFO_VAL(externalParameterIndex);
00132             MN_INFO_VAL2("internal",i);
00133 #endif
00134             good = false;
00135          }
00136       }
00137       if(!good) {
00138 #ifdef WARNINGMSG
00139          MN_ERROR_MSG("Minuit does not accept user specified Gradient. To force acceptance, override 'virtual bool CheckGradient() const' of FCNGradientBase.h in the derived class.");
00140 #endif
00141          assert(good);
00142       }
00143    }
00144    
00145    MnAlgebraicSymMatrix mat(n);
00146    double dcovar = 1.;
00147    if(st.HasCovariance()) {
00148       for(unsigned int i = 0; i < n; i++)       
00149          for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
00150       dcovar = 0.;
00151    } else {
00152       for(unsigned int i = 0; i < n; i++)       
00153          mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
00154    }
00155    MinimumError err(mat, dcovar);
00156    double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
00157    MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
00158    
00159    NegativeG2LineSearch ng2ls;
00160    if(ng2ls.HasNegativeG2(dgrad, prec)) {
00161       Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra);
00162       state = ng2ls(fcn, state, ngc, prec);
00163    }
00164    
00165    if(stra.Strategy() == 2 && !st.HasCovariance()) {
00166       //calculate full 2nd derivative
00167       MinimumState tmpState = MnHesse(stra)(fcn, state, st.Trafo());
00168       return MinimumSeed(tmpState, st.Trafo());
00169    }
00170    
00171    return MinimumSeed(state, st.Trafo());
00172 }
00173 
00174    }  // namespace Minuit2
00175 
00176 }  // namespace ROOT

Generated on Tue Jul 5 14:37:10 2011 for ROOT_528-00b_version by  doxygen 1.5.1