diff --git a/math/minuit2/inc/Minuit2/MinosError.h b/math/minuit2/inc/Minuit2/MinosError.h index 4f540c0758374..85b2ea39baad8 100644 --- a/math/minuit2/inc/Minuit2/MinosError.h +++ b/math/minuit2/inc/Minuit2/MinosError.h @@ -55,8 +55,15 @@ class MinosError { { if (AtLowerLimit()) return LowerState().Parameter(Parameter()).LowerLimit() - fMinParValue; - if (LowerValid()) - return -1. * LowerState().Error(Parameter()) * (1. + fLower.Value()); + if (LowerValid()) { + // Minos error is value - error - aopt * error where aopt is MnCross.Value() + // If value - error is below the limit, error must be truncated at limit + double err = LowerState().Error(Parameter()); + // error is truncated if over the limit + if (LowerState().Parameter(Parameter()).HasLowerLimit()) + err = std::min(err, fMinParValue - LowerState().Parameter(Parameter()).LowerLimit()); + return -1. * err * (1. + fLower.Value()); + } // return Hessian Error in case is invalid return -LowerState().Error(Parameter()); } @@ -64,8 +71,14 @@ class MinosError { { if (AtUpperLimit()) return UpperState().Parameter(Parameter()).UpperLimit() - fMinParValue; - if (UpperValid()) - return UpperState().Error(Parameter()) * (1. + fUpper.Value()); + if (UpperValid()) { + // Minos error is value + error + aopt * error where aopt is MnCross.Value() + // If value + error is over the limit, err must be truncated at limit + double err = UpperState().Error(Parameter()); + if (UpperState().Parameter(Parameter()).HasUpperLimit()) + err = std::min(err, UpperState().Parameter(Parameter()).UpperLimit() - fMinParValue); + return err * (1. + fUpper.Value()); + } // return Hessian Error in case is invalid return UpperState().Error(Parameter()); } diff --git a/math/minuit2/src/MnContours.cxx b/math/minuit2/src/MnContours.cxx index d181c52ba5ea3..95cd0956c2c25 100644 --- a/math/minuit2/src/MnContours.cxx +++ b/math/minuit2/src/MnContours.cxx @@ -40,6 +40,7 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int unsigned int nfcn = 0; MnPrint print("MnContours"); + print.Debug("MnContours: finding ",npoints," contours points for ",px,py," at level ",fFCN.Up()," from value ",fMinimum.Fval()); std::vector> result; result.reserve(npoints); @@ -47,50 +48,67 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int // double edmmax = 0.5*0.05*fFCN.Up()*1.e-3; // double toler = 0.05; - double toler = 0.1; // use same default value as in Minos - - // get first four points - // std::cout<<"MnContours: get first 4 params."< ex = meX(); + std::pair ex = mnex(); + + print.Debug("Minos error for p0: ",ex.first,ex.second); - MinosError meY = minos.Minos(py); - nfcn += meY.NFcn(); - if (!meY.IsValid()) { + MinosError mney = minos.Minos(py); + nfcn += mney.NFcn(); + if (!mney.IsValid()) { print.Error("unable to find second two points"); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } - std::pair ey = meY(); + std::pair ey = mney(); - MnMigrad migrad(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1)))); + print.Debug("Minos error for p0: ",ey.first,ey.second); - migrad.Fix(px); - migrad.SetValue(px, valx + ex.second); - FunctionMinimum exy_up = migrad(); - nfcn += exy_up.NFcn(); - if (!exy_up.IsValid()) { - print.Error("unable to find Upper y Value for x Parameter", px); - return ContoursError(px, py, result, meX, meY, nfcn); - } + // if Minos is not at limits we can use migrad to find the other corresponding point coordinate + MnMigrad migrad0(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1)))); - migrad.SetValue(px, valx + ex.first); - FunctionMinimum exy_lo = migrad(); + + // start from minimizing in p1 and fixing p0 to Minos value + migrad0.Fix(px); + migrad0.SetValue(px, valx + ex.first); + FunctionMinimum exy_lo = migrad0(); nfcn += exy_lo.NFcn(); if (!exy_lo.IsValid()) { print.Error("unable to find Lower y Value for x Parameter", px); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } + print.Debug("Minimum p1 found for p0 set to ",migrad0.Value(px)," is ",exy_lo.UserState().Value(py),"fcn = ",exy_lo.Fval()); + + migrad0.SetValue(px, valx + ex.second); + FunctionMinimum exy_up = migrad0(); + nfcn += exy_up.NFcn(); + if (!exy_up.IsValid()) { + print.Error("unable to find Upper y Value for x Parameter", px); + return ContoursError(px, py, result, mnex, mney, nfcn); + } + print.Debug("Minimum p1 found for p0 set to ",migrad0.Value(px)," is ",exy_up.UserState().Value(py),"fcn = ",exy_up.Fval()); + + MnMigrad migrad1(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1)))); migrad1.Fix(py); migrad1.SetValue(py, valy + ey.second); @@ -98,17 +116,21 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int nfcn += eyx_up.NFcn(); if (!eyx_up.IsValid()) { print.Error("unable to find Upper x Value for y Parameter", py); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } + print.Debug("Minimum p0 found for p1 set to ",migrad1.Value(py)," is ",eyx_up.UserState().Value(px),"fcn = ",eyx_up.Fval()); migrad1.SetValue(py, valy + ey.first); FunctionMinimum eyx_lo = migrad1(); nfcn += eyx_lo.NFcn(); if (!eyx_lo.IsValid()) { print.Error("unable to find Lower x Value for y Parameter", py); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } + print.Debug("Minimum p0 found for p1 set to ",migrad1.Value(py)," is ",eyx_lo.UserState().Value(px),"fcn = ",eyx_lo.Fval()); + + double scalx = 1. / (ex.second - ex.first); double scaly = 1. / (ey.second - ey.first); @@ -119,7 +141,7 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int MnUserParameterState upar = fMinimum.UserState(); - print.Info("List of found points", '\n', " Parameter x is", upar.Name(px), '\n', " Parameter y is", upar.Name(py), + print.Debug("List of first 4 found contour points", '\n', " Parameter x is", upar.Name(px), '\n', " Parameter y is", upar.Name(py), '\n', result[0], '\n', result[1], '\n', result[2], '\n', result[3]); upar.Fix(px); @@ -130,15 +152,17 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int par[1] = py; MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy); + // find the remaining points of the contour for (unsigned int i = 4; i < npoints; i++) { + // find the two neighbouring points with largest separation auto idist1 = result.end() - 1; auto idist2 = result.begin(); double dx = idist1->first - (idist2)->first; double dy = idist1->second - (idist2)->second; double bigdis = scalx * scalx * dx * dx + scaly * scaly * dy * dy; - for (auto ipair = result.begin(); ipair != result.end() - 1; ++ipair) { + for (auto ipair = result.begin(); ipair != result.end() - 1; ++ipair) { double distx = ipair->first - (ipair + 1)->first; double disty = ipair->second - (ipair + 1)->second; double dist = scalx * scalx * distx * distx + scaly * scaly * disty * disty; @@ -157,11 +181,17 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int if (nfcn > maxcalls) { print.Error("maximum number of function calls exhausted"); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } + print.Debug("Find new contour point between points with max sep: (",idist1->first,", ",idist1->second,") and (", + idist2->first,", ",idist2->second,") with weights ",a1,a2); + // find next point between the found 2 with max separation + // start from point situated at the middle (a1,a2=0.5) + // and direction double xmidcr = a1 * idist1->first + a2 * (idist2)->first; double ymidcr = a1 * idist1->second + a2 * (idist2)->second; + // direction is the perpendicular one double xdir = (idist2)->second - idist1->second; double ydir = idist1->first - (idist2)->first; double scalfac = sca * std::max(std::fabs(xdir * scalx), std::fabs(ydir * scaly)); @@ -177,30 +207,39 @@ ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int MnCross opt = cross(par, pmid, pdir, toler, maxcalls); nfcn += opt.NFcn(); if (!opt.IsValid()) { - // if(a1 > 0.5) { - if (sca < 0.) { + if(a1 > 0.5) { + // LM 20/10/23 : remove switch of direction and look instead closer (this is what is done in TMinuit) + // should we try again closer to P2 (e.g. a1=0.25, a2 = 0.75) if failing? + //if (sca < 0.) { print.Error("unable to find point on Contour", i + 1, '\n', "found only", i, "points"); - return ContoursError(px, py, result, meX, meY, nfcn); + return ContoursError(px, py, result, mnex, mney, nfcn); } - // a1 = 0.75; - // a2 = 0.25; - // std::cout<<"*****switch direction"< prec.Eps2() ? 1. / dgrad.G2()(i) : - (dgrad.G2()(i) >= 0) ? 1./prec.Eps2() : -1./prec.Eps2()); + // if G2 is small better using an arbitrary value (e.g. 1) + mat(i, i) = std::fabs(dgrad.G2()(i)) > prec.Eps() ? 1. / dgrad.G2()(i) : 1.0; } MinimumError err(mat, dcovar); @@ -177,6 +177,7 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn &fcn, const AnalyticalGradie // do this only when we have not computed the Hessian or always ? if (!computedHessian) { // check if minimum state has covariance - if not use computed G2 + // should maybe this an option, sometimes is not good to re-use existing covariance if (st.HasCovariance()) { print.Info("Using existing covariance matrix"); for (unsigned int i = 0; i < n; i++) @@ -185,10 +186,9 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn &fcn, const AnalyticalGradie dcovar = 0.; } else { for (unsigned int i = 0; i < n; i++) { - // should not use a cut-off here like 1./prec.Eps() - mat(i, i) = (std::fabs(grad.G2()(i)) > prec.Eps2() ? 1. / grad.G2()(i) - : (grad.G2()(i) >= 0) ? 1. / prec.Eps2() - : -1. / prec.Eps2()); + // if G2 is very small, better using an arbitrary value (e.g. 1.) + mat(i, i) = std::fabs(grad.G2()(i)) > prec.Eps() ? 1. / grad.G2()(i) + : 1.0; } dcovar = 1.; } diff --git a/math/minuit2/src/NegativeG2LineSearch.cxx b/math/minuit2/src/NegativeG2LineSearch.cxx index dacee38261aed..88fe8b3b6112a 100644 --- a/math/minuit2/src/NegativeG2LineSearch.cxx +++ b/math/minuit2/src/NegativeG2LineSearch.cxx @@ -121,8 +121,8 @@ MinimumState NegativeG2LineSearch::operator()(const MnFcn &fcn, const MinimumSta print.Debug("Approximate new covariance after NegativeG2LS using only G2"); MnAlgebraicSymMatrix mat(n); for (unsigned int i = 0; i < n; i++) { - mat(i, i) = (std::fabs(dgrad.G2()(i)) > prec.Eps2() ? 1. / dgrad.G2()(i) : - (dgrad.G2()(i) >= 0) ? 1./prec.Eps2() : -1./prec.Eps2()); + mat(i, i) = std::fabs(dgrad.G2()(i)) > prec.Eps() ? 1. / dgrad.G2()(i) : + 1; // use an arbitrary value (e.g. 1) } MinimumError err(mat, 1.); diff --git a/tutorials/roofit/rf601_intminuit.C b/tutorials/roofit/rf601_intminuit.C index 8bf5d2e54abc2..77b7f329d621a 100644 --- a/tutorials/roofit/rf601_intminuit.C +++ b/tutorials/roofit/rf601_intminuit.C @@ -54,9 +54,6 @@ void rf601_intminuit() // Create MINUIT interface object RooMinimizer m(*nll); - // The contour algorithm in Minuit2 does not work for this tutorial - m.setMinimizerType("Minuit"); - // Activate verbose logging of MINUIT parameter space stepping m.setVerbose(true); diff --git a/tutorials/roofit/rf601_intminuit.py b/tutorials/roofit/rf601_intminuit.py index 8cb56a0821615..169fecff430bf 100644 --- a/tutorials/roofit/rf601_intminuit.py +++ b/tutorials/roofit/rf601_intminuit.py @@ -46,9 +46,6 @@ # Create MINUIT interface object m = ROOT.RooMinimizer(nll) -# The contour algorithm in Minuit2 does not work for this tutorial -m.setMinimizerType("Minuit") - # Activate verbose logging of MINUIT parameter space stepping m.setVerbose(True)