00001
00002
00003
00004
00005
00006
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
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
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
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
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
00108 unsigned int n = st.VariableParameters();
00109 const MnMachinePrecision& prec = st.Precision();
00110
00111
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
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 }
00175
00176 }