diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index a233715e4835a..09b087a24edbe 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -103,6 +103,7 @@ The `TMVA_SOFIE_GNN` tutorials have been migrated to this workflow and produce i * The header `Rstrstream.h` is deprecated and will be removed after ROOT 6.44, use instead ``. * The headers `ZipLZMA.h`, `ZipLZ4.h` and `ZipZSTD.h` are deprecated and will be removed in ROOT 6.46, use instead the public methods in the `RZip.h` interface. * Remove outdated `TSQLFile` class, was not tested for very long time, not working with only supported sqlite DB +* The `ROOT::Minuit2::Minuit2Minimizer::SetFCN(unsigned int, std::unique_ptr)` and `GetFCN()` methods are removed. They were added for RooFit, which drives Minuit2 directly now and doesn't need them anymore. To minimize an `FCNBase` object, use the Minuit2 library interface (`MnMigrad`, `MnMinimize`, etc.) directly. The `ROOT::Minuit2::NumericalDerivator` class, which was only used by RooFit, is moved to RooFit as well and is no longer part of the public Minuit2 headers. ## Build System @@ -176,7 +177,7 @@ maps) will now obtain different, mathematically consistent values. ### Skipping identically-vanishing second derivatives in numerical Hessian evaluation -Minuit2 objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function `ROOT::Minuit2::FCNBase::SecondDerivativeAlwaysVanishes()` (with corresponding setters on `ROOT::Minuit2::FCNAdapter` and `ROOT::Minuit2::Minuit2Minimizer`). +Minuit2 objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function `ROOT::Minuit2::FCNBase::SecondDerivativeAlwaysVanishes()` (with a corresponding setter on `ROOT::Minuit2::FCNAdapter`). The numerical Hessian computation in Minuit2 (`MnHesse`) skips the finite-difference evaluations for such parameter pairs, which can speed up Hesse significantly for likelihoods with many mutually independent parameters. The parameter indices in this interface always refer to the function's own full (external) parameter space, including parameters that are fixed in the minimizer, and the advertised information must hold for all parameter values; see the `FCNBase` documentation for the full contract. diff --git a/math/mathcore/inc/Math/FitMethodFunction.h b/math/mathcore/inc/Math/FitMethodFunction.h index 63eb814ac42aa..cb6aa7526d169 100644 --- a/math/mathcore/inc/Math/FitMethodFunction.h +++ b/math/mathcore/inc/Math/FitMethodFunction.h @@ -73,12 +73,13 @@ class BasicFitMethodFunction : public FunctionType { virtual double DataElement(const double *x, unsigned int i, double *g = nullptr, double *h = nullptr, bool fullHessian = false) const = 0; // flag to indicate if full Hessian computation is supported - virtual bool HasHessian() const { return false;} + bool HasHessian() const override { return false; } /** * Computes the full Hessian. Return false if Hessian is not supported */ - virtual bool Hessian(const double * x, double * hess) const { + bool Hessian(const double *x, double *hess) const override + { //return full Hessian of the objective function which is Sum(F(i)) unsigned int np = NPoints(); unsigned int ndim = NDim(); diff --git a/math/mathcore/inc/Math/Functor.h b/math/mathcore/inc/Math/Functor.h index 8e6f848eec2c2..e7f0857ca3719 100644 --- a/math/mathcore/inc/Math/Functor.h +++ b/math/mathcore/inc/Math/Functor.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -202,10 +203,13 @@ class GradFunctor : public IGradientFunctionMultiDim { * @param f : function object computing the function value * @param dim : number of function dimension * @param g : function object computing the function gradient + * @param h : optional function object computing the function Hessian, + * filling its second argument as a row-major dim * dim array */ - GradFunctor(std::function const&f, unsigned int dim, - std::function const& g) - : fDim{dim}, fFunc{f}, fGradFunc{g} + GradFunctor(std::function const &f, unsigned int dim, + std::function const &g, + std::function const &h = nullptr) + : fDim{dim}, fFunc{f}, fGradFunc{g}, fHessFunc{h} {} // Clone of the function handler (use copy-ctor). @@ -224,6 +228,17 @@ class GradFunctor : public IGradientFunctionMultiDim { fGradFunc(x, g); } + bool HasHessian() const override { return bool(fHessFunc); } + + bool Hessian(const double *x, double *hess) const override + { + if (!fHessFunc) { + return false; + } + fHessFunc(x, hess); + return true; + } + private : inline double DoEval (const double * x) const override { @@ -246,6 +261,7 @@ private : std::function fFunc; std::function fDerivFunc; std::function fGradFunc; + std::function fHessFunc; }; @@ -285,29 +301,42 @@ class GradFunctor1D : public IGradientFunctionOneDim { the function evaluation and the other for the derivative. The member functions must take a double as argument and return a double */ - template - GradFunctor1D(const PtrObj& p, MemFn memFn, GradMemFn gradFn) + template ::value, bool> = true> + GradFunctor1D(const PtrObj &p, MemFn memFn, GradMemFn gradFn) : fFunc{std::bind(memFn, p, std::placeholders::_1)}, fDerivFunc{std::bind(gradFn, p, std::placeholders::_1)} {} - /// Specialized constructor from 2 function objects implementing double /// operator()(double x). The first one for the function evaluation and the - /// second one implementing the function derivative. - GradFunctor1D(std::function const& f, std::function const& g) - : fFunc{f}, fDerivFunc{g} + /// second one implementing the function derivative. Optionally, a third + /// function object implementing the second derivative can be passed. + GradFunctor1D(std::function const &f, std::function const &g, + std::function const &g2 = nullptr) + : fFunc{f}, fDerivFunc{g}, fSecondDerivFunc{g2} {} // clone of the function handler (use copy-ctor) GradFunctor1D * Clone() const override { return new GradFunctor1D(*this); } + bool HasHessian() const override { return bool(fSecondDerivFunc); } + private : inline double DoEval (double x) const override { return fFunc(x); } inline double DoDerivative (double x) const override { return fDerivFunc(x); } + inline double DoSecondDerivative(double x) const override + { + if (!fSecondDerivFunc) { + throw std::runtime_error( + "The second derivative evaluation is not implemented for this function (see HasHessian())"); + } + return fSecondDerivFunc(x); + } std::function fFunc; std::function fDerivFunc; + std::function fSecondDerivFunc; }; diff --git a/math/mathcore/inc/Math/IFunction.h b/math/mathcore/inc/Math/IFunction.h index bc0b029d4bceb..a24a970c130d7 100644 --- a/math/mathcore/inc/Math/IFunction.h +++ b/math/mathcore/inc/Math/IFunction.h @@ -33,6 +33,7 @@ #include "Math/IFunctionfwd.h" +#include namespace ROOT { namespace Math { @@ -111,6 +112,15 @@ namespace ROOT { Gradient(x, df); } + // Indicate whether this class supports second derivative (Hessian) calculations, + // i.e., if the Hessian() method is implemented. + virtual bool HasHessian() const { return false; } + + /// Evaluate all second derivatives (the Hessian matrix) at a point x, + /// filling hess as a row-major NDim() * NDim() array. + /// Return false if the second derivatives are not implemented (see HasHessian()). + virtual bool Hessian(const T * /*x*/, T * /*hess*/) const { return false; } + /// Return the partial derivative with respect to the passed coordinate. T Derivative(const T *x, unsigned int icoord = 0) const { return DoDerivative(x, icoord); } @@ -178,6 +188,10 @@ namespace ROOT { // if it inherits from ROOT::Math::IGradientFunctionOneDim. virtual bool HasGradient() const { return false; } + // Indicate whether this class supports second derivative calculations, + // i.e., if the SecondDerivative() method is implemented. + virtual bool HasHessian() const { return false; } + /// Return the derivative of the function at a point x /// Use the private method DoDerivative double Derivative(double x) const { return DoDerivative(x); } @@ -188,6 +202,21 @@ namespace ROOT { /// Compatibility method with multi-dimensional interface for Gradient. void Gradient(const double *x, double *g) const { g[0] = DoDerivative(*x); } + /// Return the second derivative of the function at a point x. + /// Check HasHessian() to see whether it is implemented. + double SecondDerivative(double x) const { return DoSecondDerivative(x); } + + /// Compatibility method with multi-dimensional interface for the Hessian. + /// Return false if the second derivative is not implemented (see HasHessian()). + bool Hessian(const double *x, double *hess) const + { + if (!HasHessian()) { + return false; + } + hess[0] = DoSecondDerivative(*x); + return true; + } + /// Optimized method to evaluate at the same time the function value and derivative at a point x. /// Often both value and derivatives are needed and it is often more efficient to compute them at the same time. /// Derived class should implement this method if performances play an important role and if it is faster to @@ -208,6 +237,13 @@ namespace ROOT { /// Function to evaluate the derivative with respect each coordinate. To be implemented by the derived class. virtual double DoDerivative(double) const { return 0.; } + + /// Function to evaluate the second derivative. To be implemented by the derived class (see HasHessian()). + virtual double DoSecondDerivative(double) const + { + throw std::runtime_error( + "The second derivative evaluation is not implemented for this function (see HasHessian())"); + } }; diff --git a/math/minuit2/CMakeLists.txt b/math/minuit2/CMakeLists.txt index aebf783361302..3e8f70b29406d 100644 --- a/math/minuit2/CMakeLists.txt +++ b/math/minuit2/CMakeLists.txt @@ -91,7 +91,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) Minuit2/ModularFunctionMinimizer.h Minuit2/NegativeG2LineSearch.h Minuit2/Numerical2PGradientCalculator.h - Minuit2/NumericalDerivator.h Minuit2/ParametricFunction.h Minuit2/ScanBuilder.h Minuit2/ScanMinimizer.h @@ -152,7 +151,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) src/ModularFunctionMinimizer.cxx src/NegativeG2LineSearch.cxx src/Numerical2PGradientCalculator.cxx - src/NumericalDerivator.cxx src/ParametricFunction.cxx src/ScanBuilder.cxx src/SimplexBuilder.cxx diff --git a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h index c188f6038df53..7814f14ee1ed7 100644 --- a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h +++ b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h @@ -82,15 +82,6 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer { /// set the function implementing Hessian computation void SetHessianFunction(std::function, double *)> hfunc) override; - /// Set a predicate advertising which mixed second derivatives of the - /// minimized function are identically zero, so that MnHesse can skip the - /// corresponding finite-difference evaluations. The predicate must - /// fulfill the contract documented for - /// FCNBase::SecondDerivativeAlwaysVanishes(); in particular, the indices - /// refer to the full (external) parameter space. Like SetHessianFunction, - /// this must be called after SetFunction. - void SetSecondDerivativeAlwaysVanishesFunc(std::function func); - /// set free variable bool SetVariable(unsigned int ivar, const std::string &name, double val, double step) override; @@ -288,12 +279,6 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer { /// return the minimizer state (containing values, step size , etc..) const ROOT::Minuit2::MnUserParameterState &State() { return fState; } - /// To set the function directly to a Minuit 2 function. - void SetFCN(unsigned int nDim, std::unique_ptr fcn); - - const ROOT::Minuit2::FCNBase *GetFCN() const { return fMinuitFCN.get(); } - ROOT::Minuit2::FCNBase *GetFCN() { return fMinuitFCN.get(); } - protected: // protected function for accessing the internal Minuit2 object. Needed for derived classes diff --git a/math/minuit2/src/Minuit2Minimizer.cxx b/math/minuit2/src/Minuit2Minimizer.cxx index 54b87884ae7b1..94e401325e035 100644 --- a/math/minuit2/src/Minuit2Minimizer.cxx +++ b/math/minuit2/src/Minuit2Minimizer.cxx @@ -416,15 +416,6 @@ void Minuit2Minimizer::SetHessianFunction(std::functionSetHessianFunction(hfunc); } -void Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(std::function func) -{ - // not supported for Fumili, whose FCN is not an FCNAdapter - if (fUseFumili) return; - auto fcn = static_cast(fMinuitFCN.get()); - if (!fcn) return; - fcn->SetSecondDerivativeAlwaysVanishesFunc(std::move(func)); -} - namespace { ROOT::Minuit2::MnStrategy customizedStrategy(unsigned int strategyLevel, ROOT::Math::MinimizerOptions const &options) @@ -1372,12 +1363,6 @@ void Minuit2Minimizer::SetStorageLevel(int level) fMinimizer->Builder().SetStorageLevel(level); } -void Minuit2Minimizer::SetFCN(unsigned int nDim, std::unique_ptr fcn) -{ - fDim = nDim; - fMinuitFCN = std::move(fcn); -} - } // end namespace Minuit2 } // end namespace ROOT diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index c0186c5e9bc0c..d27fe51a92221 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -475,6 +475,67 @@ void functorCodegenImpl(RooArg_t &arg, RooArgList const &variables, CodegenConte "}\n" "} // namespace clad::custom_derivatives\n"; + // If the functor also provides second derivatives, emit the + // forward-mode derivative and its pullback, which Clad needs to + // generate second derivatives of code calling the wrapper (it + // implements Hessians as forward-mode derivatives that are then + // differentiated in reverse mode). + if (arg.function()->HasHessian()) { + code += "#include \n" + "namespace clad::custom_derivatives {\n\n" + "clad::ValueAndPushforward " + + wrapperName + + "_pushforward(double const *x, double const *d_x) {\n" + " double grad[" + + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Gradient(x, grad);\n" + " double dot = 0.;\n" + " for (int i = 0; i < " + + nStr + + "; ++i) {\n" + " dot += grad[i] * d_x[i];\n" + " }\n" + " return {" + + funcAddrCasted + + "->operator()(x), dot};\n" + "}\n\n" + "void " + + wrapperName + + "_pushforward_pullback(double const *x, double const *d_x, clad::ValueAndPushforward d_y, double *d_out_x, double *d_out_d_x) {\n" + " double grad[" + + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Gradient(x, grad);\n" + " double hess[" + + nStr + " * " + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Hessian(x, hess);\n" + " for (int i = 0; i < " + + nStr + + "; ++i) {\n" + " double hdot = 0.;\n" + " for (int j = 0; j < " + + nStr + + "; ++j) {\n" + " hdot += hess[i * " + + nStr + + " + j] * d_x[j];\n" + " }\n" + " d_out_x[i] += d_y.value * grad[i] + d_y.pushforward * hdot;\n" + " d_out_d_x[i] += d_y.pushforward * grad[i];\n" + " }\n" + "}\n" + "} // namespace clad::custom_derivatives\n"; + } + gInterpreter->Declare(code.c_str()); } diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx index 29acd2c4ded64..58c1725efa9e3 100644 --- a/roofit/histfactory/test/testHistFactory.cxx +++ b/roofit/histfactory/test/testHistFactory.cxx @@ -12,17 +12,19 @@ #include #include -#include -#include -#include #include -#include -#include -#include -#include +#include +#include +#include +#include #include +#include +#include #include -#include +#include +#include +#include +#include #include #include @@ -610,8 +612,21 @@ TEST_P(HFFixtureFit, Fit) } using namespace RooFit; - std::unique_ptr fitResult{simPdf->fitTo( - *data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()), Save(), PrintLevel(verbose ? 1 : -1))}; + std::unique_ptr nll{ + simPdf->createNLL(*data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()))}; + RooMinimizer::Config cfg; + if (evalBackend == RooFit::EvalBackend::Codegen()) { + // Make sure we use both analytical gradient and Hessian + static_cast(*nll).generateGradient(); + static_cast(*nll).generateHessian(); + cfg.useGradient = true; + cfg.useHessian = true; + } + RooMinimizer minim{*nll, cfg}; + minim.setPrintLevel(verbose ? 1 : -1); + minim.minimize("Minuit2", "Migrad"); + minim.hesse(); + std::unique_ptr fitResult{minim.save()}; ASSERT_NE(fitResult, nullptr); if (verbose) fitResult->Print("v"); diff --git a/roofit/roofit/src/RooONNXFunc.cxx b/roofit/roofit/src/RooONNXFunc.cxx index b3f9e86506948..ac8a7d54b6e04 100644 --- a/roofit/roofit/src/RooONNXFunc.cxx +++ b/roofit/roofit/src/RooONNXFunc.cxx @@ -192,6 +192,7 @@ struct RooONNXFunc::RuntimeCache { RooFit::Detail::AnyWithVoidPtr _session; RooFit::Detail::AnyWithVoidPtr _d_session; + RooFit::Detail::AnyWithVoidPtr _zero_session; ///< Zero-weight session, the tangent for forward-mode derivatives. Func _func; }; @@ -295,20 +296,26 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: gInterpreter->ProcessLine(("std::size(" + namespaceName + "::inputTensorDims);").c_str())); // Per-input-tensor parameter / argument lists used by the JIT'd code below. - std::string innerParams; // "float const *input0, float const *input1, ..." - std::string innerArgs; // "input0, input1, ..." - std::string outerDoubleParams; // "double const *input0, double const *input1, ..." - std::string cladInputs; // "input0, input1, ..." (for clad::gradient param spec) + std::string innerParams; // "float const *input0, float const *input1, ..." + std::string innerArgs; // "input0, input1, ..." + std::string innerTangentParams; // "float const *d_input0, float const *d_input1, ..." + std::string innerTangentArgs; // "d_input0, d_input1, ..." + std::string outerDoubleParams; // "double const *input0, double const *input1, ..." + std::string cladInputs; // "input0, input1, ..." (for clad::gradient param spec) for (std::size_t i = 0; i < nInputTensors; ++i) { std::string istr = std::to_string(i); if (i > 0) { innerParams += ", "; innerArgs += ", "; + innerTangentParams += ", "; + innerTangentArgs += ", "; outerDoubleParams += ", "; cladInputs += ", "; } innerParams += "float const *input" + istr; innerArgs += "input" + istr; + innerTangentParams += "float const *d_input" + istr; + innerTangentArgs += "d_input" + istr; outerDoubleParams += "double const *input" + istr; cladInputs += "input" + istr; } @@ -363,6 +370,47 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: gInterpreter->ProcessLine(("clad::gradient(" + namespaceName + "::roo_wrapper, \"" + cladInputs + "\");").c_str()); + // Also generate roo_inner_wrapper_pushforward: the forward-mode derivative + // of the inner wrapper, which takes the session and input tangents as + // function arguments (requesting the derivative of the outer wrapper + // generates it, the same trick as for the pullback above). It is the basis + // of the exact second derivatives below. + gInterpreter->ProcessLine(("clad::differentiate(" + namespaceName + "::roo_wrapper, \"input0[0]\");").c_str()); + + // The session tangent for the pushforward: the weights are constants, so + // their tangent is zero (a default-constructed Session holds the actual + // weights). Intermediate tensors inside this session are used as tangent + // scratch space during the forward pass. + _runtime->_zero_session.emplace(sessionName); + auto ptrZeroSession = toPtrString(_runtime->_zero_session.ptr, sessionName); + gInterpreter->ProcessLine((ptrZeroSession + "->SetWeightsToZero();").c_str()); + + // Directional-derivative wrapper: reverse differentiation of it w.r.t. + // both the inputs and the tangent direction gives the exact + // Hessian-vector product H . d_input (in the input adjoints) and the + // gradient (in the tangent adjoints). Both must be requested as active: + // for a non-varied argument clad would look up the custom helper pullbacks + // with reduced signatures, not find them, and silently fall back to + // differentiating the BLAS calls, which yields zero. + { + std::ostringstream ss; + ss << "namespace " << namespaceName << " {\n\n" + << "float roo_dir(Session const &session, Session const &zeroSession, " << innerParams << ", " + << innerTangentParams << ") {\n" + << " return roo_inner_wrapper_pushforward(session, " << innerArgs << ", zeroSession, " << innerTangentArgs + << ").pushforward;\n" + << "}\n\n" + << "float roo_dir_outer(Session const &session, Session const &zeroSession, " << innerParams << ", " + << innerTangentParams << ") {\n" + << " return roo_dir(session, zeroSession, " << innerArgs << ", " << innerTangentArgs << ");\n" + << "}\n\n" + << "} // namespace " << namespaceName << "\n"; + gInterpreter->Declare(ss.str().c_str()); + } + gInterpreter->ProcessLine( + ("clad::gradient(" + namespaceName + "::roo_dir_outer, \"" + cladInputs + ", " + innerTangentArgs + "\");") + .c_str()); + // The codegen call site (CodegenImpl::codegenImpl(RooONNXFunc)) passes one // double-array argument per input tensor. Emit roo_outer_wrapper and the matching // custom-derivative pullback with the corresponding number of parameters. @@ -416,6 +464,108 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: << " d_input" << i << "[i] += d_inputFlt" << i << "[i];\n" << " }\n"; } + ss << "}\n\n"; + + // Custom derivatives for Clad Hessians, which are implemented in clad + // as forward-mode derivatives that are then differentiated in reverse + // mode: the forward pass needs a "pushforward" for the opaque + // roo_outer_wrapper call, and the reverse pass over it needs the + // matching "pushforward_pullback". Both are exact: they are built on + // the clad-generated roo_inner_wrapper_pushforward and roo_dir_pullback + // (see above). + + // Comma-separated parameter helper for the emitted code. + std::string tangentDoubleParams; // "double const *d_input0, ..." + for (std::size_t i = 0; i < nInputTensors; ++i) { + std::string istr = std::to_string(i); + if (i > 0) { + tangentDoubleParams += ", "; + } + tangentDoubleParams += "double const *d_input" + istr; + } + + // Emits float conversion buffers for the inputs and the tangents. + auto emitFloatBuffers = [&]() { + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " float inputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " float dInputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n" + << " inputFlt" << i << "[i] = input" << i << "[i];\n" + << " dInputFlt" << i << "[i] = d_input" << i << "[i];\n" + << " }\n"; + } + }; + + // The pushforward evaluates the function value together with the exact + // directional derivative grad . d_input, in one forward pass. + ss << "clad::ValueAndPushforward roo_outer_wrapper_pushforward(" << outerDoubleParams << ", " + << tangentDoubleParams << ") {\n" + << " using namespace ::" << namespaceName << ";\n" + << " auto &session = *" << ptrSession << ";\n" + << " auto &zeroSession = *" << ptrZeroSession << ";\n"; + emitFloatBuffers(); + ss << " auto vp = roo_inner_wrapper_pushforward(session, "; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << "inputFlt" << i << ", "; + } + ss << "zeroSession"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", dInputFlt" << i; + } + ss << ");\n" + << " return {vp.value, vp.pushforward};\n" + << "}\n\n"; + + // The pullback of the pushforward needs second derivatives only in the + // form of the Hessian-vector product H . d_input, which roo_dir_pullback + // evaluates exactly (together with the gradient, as the adjoint of the + // tangent direction). + ss << "void roo_outer_wrapper_pushforward_pullback(" << outerDoubleParams << ", " << tangentDoubleParams + << ", clad::ValueAndPushforward d_y"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", double *d_out_input" << i; + } + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", double *d_out_d_input" << i; + } + ss << ") {\n" + << " using namespace ::" << namespaceName << ";\n" + << " auto &session = *" << ptrSession << ";\n" + << " auto &zeroSession = *" << ptrZeroSession << ";\n" + << " // Fresh adjoint sessions for every call: the clad-generated\n" + << " // second-order pullback does not restore the adjoint state inside\n" + << " // the sessions to zero (clad bug, see the \"clad referenced\n" + << " // '_tracker...' before its declaration\" warnings it prints during\n" + << " // generation), so reusing them would leak state between calls.\n" + << " Session dSession1;\n" + << " Session dSession2;\n"; + emitFloatBuffers(); + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " float hvpFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n" + << " float gradFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n"; + } + ss << " roo_dir_pullback(session, zeroSession, "; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << "inputFlt" << i << ", "; + } + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << "dInputFlt" << i << ", "; + } + ss << "1.F, &dSession1, &dSession2"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", hvpFlt" << i; + } + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", gradFlt" << i; + } + ss << ");\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n" + << " d_out_input" << i << "[i] += d_y.value * gradFlt" << i << "[i] + d_y.pushforward * hvpFlt" << i + << "[i];\n" + << " d_out_d_input" << i << "[i] += d_y.pushforward * gradFlt" << i << "[i];\n" + << " }\n"; + } ss << "}\n\n" << "} // namespace " << namespaceName << "\n" << "} // namespace clad::custom_derivatives\n"; diff --git a/roofit/roofit/test/create_onnx_model.py b/roofit/roofit/test/create_onnx_model.py index b5568dc0fbce6..90b38d729bba2 100644 --- a/roofit/roofit/test/create_onnx_model.py +++ b/roofit/roofit/test/create_onnx_model.py @@ -1,19 +1,18 @@ import numpy as np import torch import torch.nn as nn -import torch.optim as optim class SmallMLP(nn.Module): """Single-input MLP regression model.""" - def __init__(self, in_features=10, hidden=32, out_features=1): + def __init__(self, in_features=10, hidden=32, out_features=1, activation=nn.ReLU): super().__init__() self.net = nn.Sequential( nn.Linear(in_features, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, out_features), ) @@ -24,13 +23,13 @@ def forward(self, x): class TwoInputMLP(nn.Module): """Dual-input MLP: concatenates two input tensors then runs an MLP.""" - def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1): + def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.ReLU): super().__init__() self.net = nn.Sequential( nn.Linear(in_features_a + in_features_b, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, out_features), ) @@ -68,7 +67,7 @@ def export_model( return model -def run_inference_and_save(model, example_inputs, name): +def run_inference_and_save(model, example_inputs, name, with_hessian=False): """Run forward+backward with fixed inputs and save prediction + gradients.""" model = model.cpu() inputs = [t.clone().detach().requires_grad_(True) for t in example_inputs] @@ -83,6 +82,19 @@ def run_inference_and_save(model, example_inputs, name): for i, x in enumerate(inputs): np.savetxt(f"{name}_grad_{i}.txt", x.grad.detach().numpy()) + if with_hessian: + # Full Hessian over the concatenation of all (flattened) input + # tensors, matching the parameter ordering used on the RooFit side. + blocks = torch.autograd.functional.hessian( + lambda *args: model(*args).squeeze(), tuple(t.clone().detach() for t in example_inputs) + ) + n_per_input = [t.numel() for t in example_inputs] + hess = np.block([ + [blocks[i][j].reshape(n_per_input[i], n_per_input[j]).numpy() for j in range(len(example_inputs))] + for i in range(len(example_inputs)) + ]) + np.savetxt(f"{name}_hessian.txt", hess) + def main(): # ---- Model 1: single input ---- @@ -114,6 +126,39 @@ def main(): name="regression_mlp_two_input", ) + # ---- Models 3 and 4: like models 1 and 2, but with a smooth activation + # function so that the Hessian with respect to the inputs is non-trivial + # (a ReLU MLP is piecewise-linear in its inputs, so its input Hessian + # vanishes almost everywhere). Used for the Hessian tests. + model3 = SmallMLP(in_features=10, hidden=32, out_features=1, activation=nn.Tanh) + model3 = export_model( + model3, + input_shapes=[(10,)], + onnx_path="regression_mlp_tanh.onnx", + ) + run_inference_and_save( + model3, + example_inputs=[torch.tensor([[0.1] * 10])], + name="regression_mlp_tanh", + with_hessian=True, + ) + + model4 = TwoInputMLP(in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.Tanh) + model4 = export_model( + model4, + input_shapes=[(10,), (5,)], + onnx_path="regression_mlp_two_input_tanh.onnx", + ) + run_inference_and_save( + model4, + example_inputs=[ + torch.tensor([[0.1] * 10]), + torch.tensor([[0.2] * 5]), + ], + name="regression_mlp_two_input_tanh", + with_hessian=True, + ) + if __name__ == "__main__": main() diff --git a/roofit/roofit/test/testRooONNXFunc.cxx b/roofit/roofit/test/testRooONNXFunc.cxx index 06b2638a84ffe..f2f2d5975261e 100644 --- a/roofit/roofit/test/testRooONNXFunc.cxx +++ b/roofit/roofit/test/testRooONNXFunc.cxx @@ -2,6 +2,7 @@ // Authors: Jonas Rembser, CERN 2026 #include +#include #include #include #include @@ -197,4 +198,119 @@ TEST(RooONNXFunc, Basic_CodegenAD_2Tensors) EXPECT_NEAR(output_vec[10 + i], refGrad1[i], 1e-5); } } + +/// Validate the Clad Hessian of a RooONNXFunc against a PyTorch reference. +/// The Hessian tests use models with tanh activations: a ReLU MLP is +/// piecewise-linear in its inputs, so its input Hessian vanishes almost +/// everywhere and would not validate anything. +TEST(RooONNXFunc, CodegenHessian) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0]; + std::vector refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt"); + std::vector refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt"); + + const std::size_t n = 10; + + RooArgList args; + fillArgs(args, n); + + RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5); + + roo_final.generateGradient(); + + std::vector grad(n); + roo_final.gradient(grad.data()); + for (std::size_t i = 0; i < n; ++i) { + EXPECT_NEAR(grad[i], refGrad[i], 1e-5); + } + + roo_final.generateHessian(); + + // The second derivatives are exact Hessian-vector products (reverse-mode + // derivatives of the forward-mode derivative of the generated code), so + // they merit the same tolerance as the gradient checks. + std::vector hess(n * n); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < n * n; ++i) { + EXPECT_NEAR(hess[i], refHess[i], 1e-5); + } +} + +/// Test the Clad Hessian of a RooONNXFunc with two input tensors, including +/// the Hessian blocks that mix the two tensors. +TEST(RooONNXFunc, CodegenHessian_2Tensors) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_two_input_tanh_pred.txt")[0]; + std::vector refHess = readDoublesFromFile("regression_mlp_two_input_tanh_hessian.txt"); + + RooArgList args0; + fillArgs(args0, 10, 0.1, "a"); + RooArgList args1; + fillArgs(args1, 5, 0.2, "b"); + + RooONNXFunc roo_func{"func", "", {args0, args1}, "regression_mlp_two_input_tanh.onnx"}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5); + + roo_final.generateHessian(); + + const std::size_t nTotal = 10 + 5; + std::vector hess(nTotal * nTotal); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < nTotal * nTotal; ++i) { + EXPECT_NEAR(hess[i], refHess[i], 1e-5); + } +} + +/// Test the Clad Hessian of a compound expression: the square of a +/// RooONNXFunc. When the RooONNXFunc is the top-level function, the adjoint +/// of the primal value in the emitted pushforward pullback stays zero; the +/// product rule in H(f^2) = 2 * (f * H + grad * grad^T) exercises it. +TEST(RooONNXFunc, CodegenHessianCompound) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0]; + std::vector refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt"); + std::vector refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt"); + + const std::size_t n = 10; + + RooArgList args; + fillArgs(args, n); + + RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"}; + RooProduct square{"square", "", {roo_func, roo_func}}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{square, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred * refPred, 1e-5); + + roo_final.generateHessian(); + + std::vector hess(n * n); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < n; ++i) { + for (std::size_t j = 0; j < n; ++j) { + const double ref = 2. * (refPred * refHess[i * n + j] + refGrad[i] * refGrad[j]); + EXPECT_NEAR(hess[i * n + j], ref, 1e-5); + } + } +} #endif diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 4e2fe10fcd4b5..3cef3d2e82b80 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -14,12 +14,14 @@ if(NOT WIN32) src/TestStatistics/LikelihoodGradientJob.cxx src/TestStatistics/LikelihoodJob.cxx src/TestStatistics/MinuitFcnGrad.cxx + src/TestStatistics/NumericalDerivator.cxx ) set(RooFitMPTestStatisticsHeaders src/TestStatistics/LikelihoodGradientJob.h src/TestStatistics/LikelihoodSerial.h src/TestStatistics/MinuitFcnGrad.h src/TestStatistics/LikelihoodJob.h + src/TestStatistics/NumericalDerivator.h ) list(APPEND EXTRA_LIBRARIES RooFitMultiProcess) #FIXME: The ProcessTimer.h exposes json in its interface: @@ -486,8 +488,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore ${EXTRA_DICT_OPTS} ) -# Minuit2 is needed in the implementation of RooMinimizerFcn, to advertise -# vanishing second derivatives to Minuit2Minimizer. +# RooMinimizer drives Minuit2 directly via its library interface, without +# going through the ROOT::Math::Minimizer plugin layer. target_link_libraries(RooFitCore PRIVATE Minuit2) # The following definitions are PUBLIC so they can also be used in ROOT-internal tests diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 4fae9157f263d..15ef5b00af88e 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -49,7 +49,6 @@ class RooMinimizer : public TObject { struct FitResult { FitResult() = default; - FitResult(const ROOT::Fit::FitConfig &fconfig); double error(unsigned int i) const { return (i < fErrors.size()) ? fErrors[i] : 0; } double lowerError(unsigned int i) const; @@ -58,6 +57,7 @@ class RooMinimizer : public TObject { double Edm() const { return fEdm; } bool IsValid() const { return fValid; } int Status() const { return fStatus; } + int CovMatrixStatus() const { return fCovStatus; } void GetCovarianceMatrix(TMatrixDSym &cov) const; bool isParameterFixed(unsigned int ipar) const; @@ -120,6 +120,9 @@ class RooMinimizer : public TObject { } ROOT::Fit::FitConfig &Config() const { return *_config; } + /// The underlying ROOT::Math::Minimizer. Returns `nullptr` when the + /// minimizer type is Minuit2, because RooMinimizer then uses Minuit2 + /// directly without going through the ROOT::Math::Minimizer interface. ROOT::Math::Minimizer *GetMinimizer() const { return _minimizer; } const FitResult &Result() const { return *_result; } @@ -183,7 +186,8 @@ class RooMinimizer : public TObject { int evalCounter() const; void zeroEvalCount(); - /// Return underlying ROOT fitter object + /// Return an object giving access to the fit configuration and the + /// internal fit result, for backwards compatibility. inline auto fitter() { return std::make_unique(&_config, _minimizer.get(), _result.get()); } int getNPar() const; @@ -215,6 +219,7 @@ class RooMinimizer : public TObject { int exec(std::string const &algoName, std::string const &statusName); bool fitFCN(); + bool runMinimizer(); bool calculateHessErrors(); bool calculateMinosErrors(); @@ -229,10 +234,25 @@ class RooMinimizer : public TObject { void fillCorrMatrix(RooFitResult &fitRes); void updateErrors(); + // Direct interface to Minuit2 that bypasses ROOT::Math::Minimizer. It is + // used whenever the minimizer type is "Minuit2". + struct Minuit2State; + bool useMinuit2Directly() const { return _cfg.minimizerType == "Minuit2"; } + void initMinuit2(); + bool minuit2Minimize(); + bool minuit2Hesse(); + bool minuit2Minos(unsigned int index, double &errLow, double &errUp); + int runMinuit2Minos(unsigned int index, double &errLow, double &errUp); + bool + minuit2Contour(unsigned int ipar, unsigned int jpar, unsigned int npoints, double errorDef, double *x, double *y); + void setMinuit2ErrorDef(double up); + void fillResultFromMinuit2(bool isValid); + RooAbsReal &_function; ROOT::Fit::FitConfig _config; ///< fitter configuration (options and parameter settings) std::unique_ptr _result; /// _minimizer; /// _minimizer; /// _minuit2; /// & /// Transfer MINUIT fit results back into RooFit objects. void RooAbsMinimizerFcn::BackProp() { - auto const &results = _context->fitter()->Result(); + auto const &results = *_context->_result; for (std::size_t index = 0; index < getNDim(); index++) { diff --git a/roofit/roofitcore/src/RooAbsMinimizerFcn.h b/roofit/roofitcore/src/RooAbsMinimizerFcn.h index 2607f9e0bc20b..cb91ad7196e0c 100644 --- a/roofit/roofitcore/src/RooAbsMinimizerFcn.h +++ b/roofit/roofitcore/src/RooAbsMinimizerFcn.h @@ -17,29 +17,35 @@ #ifndef ROO_ABS_MINIMIZER_FCN #define ROO_ABS_MINIMIZER_FCN -#include "Math/IFunction.h" - -#include "TMatrixDSym.h" - #include "RooAbsReal.h" #include "RooArgList.h" #include "RooMinimizer.h" #include "RooRealVar.h" -#include +#include +#include +#include #include #include #include #include // unique_ptr -class RooAbsMinimizerFcn { +/// Base class for the functions that RooMinimizer hands to Minuit. It +/// implements the ROOT::Minuit2::FCNBase interface, so that RooMinimizer can +/// drive Minuit2 directly. The function evaluation itself is implemented by +/// the derived classes. +class RooAbsMinimizerFcn : public ROOT::Minuit2::FCNBase { public: RooAbsMinimizerFcn(RooArgList paramList, RooMinimizer *context); - virtual ~RooAbsMinimizerFcn() = default; + ~RooAbsMinimizerFcn() override = default; - virtual void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) = 0; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double Up() const override { return _errorDef; } + void SetErrorDef(double up) override { _errorDef = up; } + /// @} /// Informs Minuit through its parameter_settings vector of RooFit parameter properties. bool synchronizeParameterSettings(std::vector ¶meters); @@ -83,6 +89,7 @@ class RooAbsMinimizerFcn { virtual void setOffsetting(bool flag) = 0; RooMinimizer::Config const &cfg() const { return _context->_cfg; } + ROOT::Math::MinimizerOptions const &minimizerOptions() const { return _context->_config.MinimizerOptions(); } inline RooRealVar &floatableParam(std::size_t i) const { @@ -122,6 +129,9 @@ class RooAbsMinimizerFcn { std::vector _floatableParamIndices; std::ofstream *_logfile = nullptr; + +private: + double _errorDef = 1.; ///< Minuit error definition, see FCNBase::Up(). }; #endif diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 133b70f7ffc77..0df6fff8999e1 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -18,10 +18,11 @@ \class RooMinimizer \ingroup Roofitcore -Wrapper class around ROOT::Math::Minimizer that -provides a seamless interface between the minimizer functionality +RooMinimizer provides a seamless interface between the minimizer functionality and the native RooFit interface. -By default the Minimizer is Minuit 2. +By default the Minimizer is Minuit 2, which RooMinimizer drives directly via the +Minuit2 library interface. Other minimizer types are reached via the +ROOT::Math::Minimizer plugin interface. RooMinimizer can minimize any RooAbsReal function with respect to its parameters. Usual choices for minimization are the object returned by RooAbsPdf::createNLL() or RooAbsReal::createChi2(). @@ -64,11 +65,26 @@ Various methods are available to control verbosity or profiling. #include "RooFitImplHelpers.h" #include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include +#include +#include #include #include #include // logic_error @@ -138,8 +154,132 @@ std::unique_ptr setOperModesDirty(RooAbsReal &function) return {}; } +// Set the global Minuit2 print level for the duration of a Minuit2 call. +class Minuit2PrintLevelRAII { +public: + Minuit2PrintLevelRAII(int level) : _prevLevel{ROOT::Minuit2::MnPrint::SetGlobalLevel(level)} {} + ~Minuit2PrintLevelRAII() { ROOT::Minuit2::MnPrint::SetGlobalLevel(_prevLevel); } + +private: + int _prevLevel; +}; + +// The Minuit2-specific extra options, either set on the minimizer options +// explicitly or taken from the global defaults for Minuit2. +ROOT::Math::IOptions const *minuit2ExtraOptions(ROOT::Math::MinimizerOptions const &options) +{ + if (auto *opts = options.ExtraOptions()) { + return opts; + } + return ROOT::Math::MinimizerOptions::FindDefault("Minuit2"); +} + +// Build the Minuit2 strategy from the strategy level and the Minuit2-specific +// extra options. +ROOT::Minuit2::MnStrategy makeMinuit2Strategy(ROOT::Math::MinimizerOptions const &options) +{ + ROOT::Minuit2::MnStrategy st{static_cast(options.Strategy())}; + ROOT::Math::IOptions const *minuit2Opt = minuit2ExtraOptions(options); + if (!minuit2Opt) { + return st; + } + auto customize = [&minuit2Opt](const char *name, auto val) { + minuit2Opt->GetValue(name, val); + return val; + }; + st.SetGradientNCycles(customize("GradientNCycles", int(st.GradientNCycles()))); + st.SetHessianNCycles(customize("HessianNCycles", int(st.HessianNCycles()))); + st.SetHessianGradientNCycles(customize("HessianGradientNCycles", int(st.HessianGradientNCycles()))); + + st.SetGradientTolerance(customize("GradientTolerance", st.GradientTolerance())); + st.SetGradientStepTolerance(customize("GradientStepTolerance", st.GradientStepTolerance())); + st.SetHessianStepTolerance(customize("HessianStepTolerance", st.HessianStepTolerance())); + st.SetHessianG2Tolerance(customize("HessianG2Tolerance", st.HessianG2Tolerance())); + + st.SetHessianCentralFDMixedDerivatives( + customize("HessianCentralFDMixedDerivatives", int(st.HessianCentralFDMixedDerivatives()))); + st.SetHessianForcePosDef(customize("HessianForcePosDef", int(st.HessianForcePosDef()))); + + return st; +} + +// Create the Minuit2 minimization algorithm for the given algorithm name. +// Returns nullptr for algorithms that RooFit can't use. +std::unique_ptr makeMinuit2Minimizer(std::string algo) +{ + std::transform(algo.begin(), algo.end(), algo.begin(), [](unsigned char c) { return std::tolower(c); }); + using namespace ROOT::Minuit2; + if (algo == "simplex") + return std::make_unique(); + if (algo == "minimize") + return std::make_unique(); + if (algo == "scan") + return std::make_unique(); + if (algo == "bfgs") + return std::make_unique(VariableMetricMinimizer::BFGSType()); + if (algo == "fumili" || algo == "fumili2") { + // Fumili needs a FumiliFCNBase, which RooFit does not provide. + return nullptr; + } + // MIGRAD is the default algorithm, also for unknown algorithm names. + return std::make_unique(); +} + +bool isFixedOrConst(ROOT::Minuit2::MinuitParameter const &par) +{ + return par.IsFixed() || par.IsConst(); +} + +// Covariance matrix element in the external parameter space, zero for fixed +// parameters and when no covariance matrix is available. +double minuit2Covariance(ROOT::Minuit2::MnUserParameterState const &state, unsigned int i, unsigned int j) +{ + if (!state.HasCovariance() || isFixedOrConst(state.Parameter(i)) || isFixedOrConst(state.Parameter(j))) { + return 0.; + } + return state.Covariance()(state.IntOfExt(i), state.IntOfExt(j)); +} + +// Covariance matrix status code with the convention of ROOT::Math::Minimizer: +// -1: not available (inversion failed or Hesse failed) +// 0: available but not positive defined +// 1: covariance only approximate +// 2: full matrix but forced positive definite +// 3: full accurate matrix +int minuit2CovMatrixStatus(ROOT::Minuit2::FunctionMinimum const *minimum, + ROOT::Minuit2::MnUserParameterState const &state) +{ + if (!minimum) { + return state.CovarianceStatus(); + } + if (minimum->HasAccurateCovar()) + return 3; + if (minimum->HasMadePosDefCovar()) + return 2; + if (minimum->HasValidCovariance()) + return 1; + if (minimum->HasCovariance()) + return 0; + return -1; +} + } // namespace +/// State of the direct Minuit2 interface. +struct RooMinimizer::Minuit2State { + /// Minuit2 parameter state: the starting point for the next minimization, + /// and the result of the last operation. + ROOT::Minuit2::MnUserParameterState state; + /// Result of the last minimization, needed for Hesse, Minos and contours. + std::unique_ptr minimum; + /// Status code with the same convention as ROOT::Minuit2::Minuit2Minimizer. + int status = 0; + /// Status of the last Minos run. + int minosStatus = -1; + /// Whether the parameter state was initialized from the parameter settings. + bool initialized = false; +}; + //////////////////////////////////////////////////////////////////////////////// /// Construct MINUIT interface to given function. Function can be anything, /// but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2 @@ -211,6 +351,7 @@ void RooMinimizer::initMinimizerFirstPart() { RooSentinel::activate(); setMinimizerType(""); + _minuit2 = std::make_unique(); _config.SetMinimizer(_cfg.minimizerType.c_str()); setEps(1.0); // default tolerance @@ -428,7 +569,7 @@ int RooMinimizer::exec(std::string const &algoName, std::string const &statusNam int RooMinimizer::hesse() { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::hesse: Error, run Migrad before Hesse!" << std::endl; _status = -1; return _status; @@ -445,7 +586,7 @@ int RooMinimizer::hesse() int RooMinimizer::minos() { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!" << std::endl; _status = -1; return _status; @@ -462,7 +603,7 @@ int RooMinimizer::minos() int RooMinimizer::minos(const RooArgSet &minosParamList) { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!" << std::endl; _status = -1; } else if (!minosParamList.empty()) { @@ -570,7 +711,7 @@ void RooMinimizer::optimizeConst(int /*flag*/) {} RooFit::OwningPtr RooMinimizer::save(const char *userName, const char *userTitle) { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::save: Error, run minimization before!" << std::endl; return nullptr; } @@ -584,7 +725,7 @@ RooFit::OwningPtr RooMinimizer::save(const char *userName, const c fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()); fitRes->setStatus(_status); - fitRes->setCovQual(_minimizer->CovMatrixStatus()); + fitRes->setCovQual(_result->fCovStatus); fitRes->setMinNLL(_result->fVal - _fcn->getOffset()); fitRes->setEDM(_result->fEdm); @@ -628,7 +769,7 @@ void RooMinimizer::fillCorrMatrix(RooFitResult &fitRes) const std::size_t nParams = _fcn->getNDim(); TMatrixDSym corrs(nParams); TMatrixDSym covs(nParams); - std::vector globalCC = _minimizer->GlobalCC(); + std::vector globalCC = _result->fGlobalCC; globalCC.resize(nParams); // pad with zeros for (std::size_t ic = 0; ic < nParams; ic++) { for (std::size_t ii = 0; ii < nParams; ii++) { @@ -684,15 +825,23 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do TMarker *point = new TMarker(var1.getVal(), var2.getVal(), 8); frame->addObject(point); - // check first if a inimizer is available. If not means - // the minimization is not done , so do it - if (_minimizer == nullptr) { + // check first if a minimization was done + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::contour: Error, run Migrad before contours!" << std::endl; return frame; } // remember our original value of ERRDEF - double errdef = _minimizer->ErrorDef(); + const double errdef = _config.MinimizerOptions().ErrorDef(); + + // compute a contour at the given error level + auto computeContour = [&](double up, double *xcoor, double *ycoor) { + if (useMinuit2Directly()) { + return minuit2Contour(index1, index2, npoints, up, xcoor, ycoor); + } + _minimizer->SetErrorDef(up); + return _minimizer->Contour(index1, index2, npoints, xcoor, ycoor); + }; double n[6]; n[0] = n1; @@ -706,13 +855,10 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do for (int ic = 0; ic < 6; ic++) { if (n[ic] > 0) { - // set the value corresponding to an n1-sigma contour - _minimizer->SetErrorDef(n[ic] * n[ic] * errdef); - - // calculate and draw the contour + // calculate and draw the contour at the level corresponding to n-sigma std::vector xcoor(npoints + 1); std::vector ycoor(npoints + 1); - bool ret = _minimizer->Contour(index1, index2, npoints, xcoor.data(), ycoor.data()); + bool ret = computeContour(n[ic] * n[ic] * errdef, xcoor.data(), ycoor.data()); if (!ret) { coutE(Minimization) << "RooMinimizer::contour(" << GetName() @@ -734,7 +880,11 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do } // restore the original ERRDEF - _minimizer->SetErrorDef(errdef); + if (useMinuit2Directly()) { + setMinuit2ErrorDef(errdef); + } else { + _minimizer->SetErrorDef(errdef); + } // restore parameter values params.assign(paramSave); @@ -805,7 +955,7 @@ RooFit::OwningPtr RooMinimizer::lastMinuitFit() // Import the results of the last fit performed, interpreting // the fit parameters as the given varList of parameters. - if (_minimizer == nullptr) { + if (_result == nullptr) { oocoutE(nullptr, InputArguments) << "RooMinimizer::save: Error, run minimization before!" << std::endl; return nullptr; } @@ -851,7 +1001,7 @@ RooFit::OwningPtr RooMinimizer::lastMinuitFit() res->setFinalParList(floatPars); res->setMinNLL(_result->fVal); res->setEDM(_result->fEdm); - res->setCovQual(_minimizer->CovMatrixStatus()); + res->setCovQual(_result->fCovStatus); res->setStatus(_result->fStatus); fillCorrMatrix(*res); @@ -954,10 +1104,7 @@ bool RooMinimizer::fitFCN() if (nPdfs == 0) { coutI(Minimization) << "[fitFCN] No discrete parameters, performing continuous minimization only" << std::endl; FreezeDisconnectedParametersRAII freeze(this, *_fcn); - bool isValid = _minimizer->Minimize(); - if (!_result) - _result = std::make_unique(); - fillResult(isValid); + bool isValid = runMinimizer(); if (isValid) updateFitConfig(); return isValid; @@ -994,12 +1141,12 @@ bool RooMinimizer::fitFCN() pdfIndices[i]->setConstant(true); } FreezeDisconnectedParametersRAII freeze(this, *_fcn); - _minimizer->Minimize(); + runMinimizer(); for (size_t i = 0; i < nPdfs; ++i) pdfIndices[i]->setConstant(wasConst[i]); - double val = _minimizer->MinValue(); + double val = _result->fVal; tried.insert(combo); nllMap[combo] = val; @@ -1016,7 +1163,7 @@ bool RooMinimizer::fitFCN() } FreezeDisconnectedParametersRAII freeze(this, *_fcn); - _minimizer->Minimize(); + runMinimizer(); coutI(Minimization) << "All NLL Values per Combination:" << std::endl; for (const auto &entry : nllMap) { @@ -1046,13 +1193,24 @@ bool RooMinimizer::fitFCN() coutI(Minimization) << ssBest.str() << std::endl; - if (!_result) - _result = std::make_unique(); - fillResult(true); + _result->fValid = true; updateFitConfig(); return true; } + +/// Run the minimization algorithm from the current minimizer state and store +/// the outcome in the internal fit result. Returns whether a valid minimum was +/// found. +bool RooMinimizer::runMinimizer() +{ + bool isValid = useMinuit2Directly() ? minuit2Minimize() : _minimizer->Minimize(); + if (!_result) + _result = std::make_unique(); + fillResult(isValid); + return isValid; +} + bool RooMinimizer::calculateHessErrors() { // compute the Hesse errors according to configuration @@ -1060,22 +1218,20 @@ bool RooMinimizer::calculateHessErrors() auto operModeRAII = setOperModesDirty(_function); - // update minimizer (recreate if not done or if name has changed - if (!updateMinimizerOptions()) { - coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; - return false; + bool ret = false; + if (useMinuit2Directly()) { + ret = minuit2Hesse(); + } else { + // update minimizer (recreate if not done or if name has changed + if (!updateMinimizerOptions()) { + coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; + return false; + } + ret = _minimizer->Hesse(); } - - // run Hesse - bool ret = _minimizer->Hesse(); if (!ret) coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error when calculating Hessian" << std::endl; - // update minimizer results with what comes out from Hesse - // in case is empty - create from a FitConfig - if (_result->fParams.empty()) - _result = std::make_unique(_config); - // re-give a minimizer instance in case it has been changed ret |= update(ret); @@ -1097,8 +1253,8 @@ bool RooMinimizer::calculateMinosErrors() auto operModeRAII = setOperModesDirty(_function); // update minimizer (but cannot re-create in this case). Must use an existing one - if (!updateMinimizerOptions(false)) { - coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; + if (!useMinuit2Directly() && !updateMinimizerOptions(false)) { + coutE(Minimization) << "RooMinimizer::calculateMinosErrors() Error re-initializing the minimizer" << std::endl; return false; } @@ -1119,9 +1275,17 @@ bool RooMinimizer::calculateMinosErrors() for (int i = 0; i < iparMax; ++i) { double elow, eup; unsigned int index = (!ipars.empty()) ? ipars[i] : i; - bool ret = _minimizer->GetMinosError(index, elow, eup); + bool ret = false; + int minosStatus = 0; + if (useMinuit2Directly()) { + ret = minuit2Minos(index, elow, eup); + minosStatus = _minuit2->minosStatus; + } else { + ret = _minimizer->GetMinosError(index, elow, eup); + minosStatus = _minimizer->MinosStatus(); + } // flags case when a new minimum has been found - if ((_minimizer->MinosStatus() & 8) != 0) { + if ((minosStatus & 8) != 0) { iparNewMin = i; } if (ret) @@ -1147,8 +1311,20 @@ bool RooMinimizer::calculateMinosErrors() void RooMinimizer::initMinimizer() { + if (useMinuit2Directly()) { + initMinuit2(); + return; + } + // Only RooMinimizerFcn supports minimizers other than Minuit2: the + // multiprocess MinuitFcnGrad is tied to Minuit2, which the constructor + // enforces. + auto *fcn = dynamic_cast(_fcn.get()); + if (!fcn) { + throw std::logic_error("RooMinimizer: minimizer type " + _cfg.minimizerType + + " is only supported with the default (non-parallel) likelihood function"); + } _minimizer = std::unique_ptr(_config.CreateMinimizer()); - _fcn->initMinimizer(*_minimizer, this); + fcn->initMinimizer(*_minimizer); _minimizer->SetVariables(_config.ParamsSettings().begin(), _config.ParamsSettings().end()); if (_cfg.setInitialCovariance) { @@ -1203,36 +1379,12 @@ void RooMinimizer::updateFitConfig() } } -RooMinimizer::FitResult::FitResult(const ROOT::Fit::FitConfig &fconfig) - : fStatus(-99), // use this special convention to flag it when printing result - fCovStatus(0), - fParams(fconfig.NPar()), - fErrors(fconfig.NPar()) -{ - // create a Fit result from a fit config (i.e. with initial parameter values - // and errors equal to step values - // The model function is NULL in this case - - // set minimizer type and algorithm - fMinimType = fconfig.MinimizerType(); - // append algorithm name for minimizer that support it - if ((fMinimType.find("Fumili") == std::string::npos) && (fMinimType.find("GSLMultiFit") == std::string::npos)) { - if (!fconfig.MinimizerAlgoType().empty()) - fMinimType += " / " + fconfig.MinimizerAlgoType(); - } - - // get parameter values and errors (step sizes) - for (unsigned int i = 0; i < fconfig.NPar(); ++i) { - const ROOT::Fit::ParameterSettings &par = fconfig.ParSettings(i); - fParams[i] = par.Value(); - fErrors[i] = par.StepSize(); - if (par.IsFixed()) - fFixedParams[i] = true; - } -} - void RooMinimizer::fillResult(bool isValid) { + if (useMinuit2Directly()) { + fillResultFromMinuit2(isValid); + return; + } ROOT::Math::Minimizer &min = *_minimizer; ROOT::Fit::FitConfig const &fconfig = _config; @@ -1270,6 +1422,7 @@ void RooMinimizer::fillResult(bool isValid) // if minimizer provides error provides also error matrix // clear in case of re-filling an existing result _result->fCovMatrix.clear(); + _result->fGlobalCC.clear(); if (min.Errors() != nullptr) { updateErrors(); @@ -1278,6 +1431,10 @@ void RooMinimizer::fillResult(bool isValid) bool RooMinimizer::update(bool isValid) { + if (useMinuit2Directly()) { + fillResultFromMinuit2(isValid); + return true; + } ROOT::Math::Minimizer &min = *_minimizer; ROOT::Fit::FitConfig const &fconfig = _config; @@ -1312,6 +1469,7 @@ void RooMinimizer::updateErrors() _result->fErrors.resize(npar); std::copy(min.Errors(), min.Errors() + npar, _result->fErrors.begin()); + _result->fGlobalCC = min.GlobalCC(); if (_result->fCovStatus != 0) { @@ -1327,6 +1485,574 @@ void RooMinimizer::updateErrors() // minos errors are set separately when calling Fitter::CalculateMinosErrors() } +//////////////////////////////////////////////////////////////////////////////// +/// Initialize the Minuit2 parameter state from the parameter settings, and +/// forget about any previous function minimum. + +void RooMinimizer::initMinuit2() +{ + using ROOT::Minuit2::MnUserParameterState; + + _minuit2->state = MnUserParameterState{}; + _minuit2->minimum.reset(); + _minuit2->status = 0; + _minuit2->minosStatus = -1; + + MnUserParameterState &state = _minuit2->state; + + for (ROOT::Fit::ParameterSettings const &par : _config.ParamsSettings()) { + if (par.IsFixed()) { + // Fixed parameters still need a step size, otherwise Minuit2 would + // treat them as constants that can't be released anymore. + const double step = par.Value() != 0 ? 0.1 * std::abs(par.Value()) : 0.1; + state.Add(par.Name(), par.Value(), step); + state.Fix(par.Name()); + continue; + } + if (par.StepSize() <= 0) { + // A parameter without a valid step size is treated as a constant. + state.Add(par.Name(), par.Value()); + } else { + state.Add(par.Name(), par.Value(), par.StepSize()); + } + if (par.IsDoubleBound()) { + state.SetLimits(par.Name(), par.LowerLimit(), par.UpperLimit()); + } else if (par.HasLowerLimit()) { + state.SetLowerLimit(par.Name(), par.LowerLimit()); + } else if (par.HasUpperLimit()) { + state.SetUpperLimit(par.Name(), par.UpperLimit()); + } + } + + if (_cfg.setInitialCovariance) { + // Diagonal covariance matrix from the parameter errors, in the packed + // lower-triangular storage of MnUserCovariance. + const unsigned int n = _fcn->getNDim(); + ROOT::Minuit2::MnUserCovariance cov{n}; + for (unsigned int i = 0; i < n; ++i) { + const double err = _fcn->floatableParam(i).getError(); + cov(i, i) = err * err; + } + state.AddCovariance(cov); + } + + _minuit2->initialized = true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Set the error definition on the function and on the function minimum, if +/// there is one. + +void RooMinimizer::setMinuit2ErrorDef(double up) +{ + _fcn->SetErrorDef(up); + if (_minuit2->minimum && _minuit2->minimum->Up() != up) { + _minuit2->minimum->SetErrorDef(up); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run the configured Minuit2 minimization algorithm, starting from the current +/// parameter state. Returns whether a valid minimum was found. + +bool RooMinimizer::minuit2Minimize() +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + + if (!_minuit2->initialized) { + initMinuit2(); + } + + auto minimizer = makeMinuit2Minimizer(options.MinimizerAlgorithm()); + if (!minimizer) { + coutE(Minimization) << "RooMinimizer::minimize: the Minuit2 algorithm \"" << options.MinimizerAlgorithm() + << "\" is not supported by RooFit" << std::endl; + return false; + } + + MnUserParameterState &state = _minuit2->state; + + // delete result of previous minimization + _minuit2->minimum.reset(); + + const unsigned int maxfcn = options.MaxFunctionCalls(); + const double tol = options.Tolerance(); + const int printLevel = options.PrintLevel(); + _fcn->SetErrorDef(options.ErrorDef()); + + if (printLevel >= 1) { + // print the real number of maxfcn used (defined in ModularFunctionMinimizer) + unsigned int maxfcnUsed = maxfcn; + if (maxfcnUsed == 0) { + const unsigned int nvar = state.VariableParameters(); + maxfcnUsed = 200 + 100 * nvar + 5 * nvar * nvar; + } + std::cout << "RooMinimizer: Minuit2 minimize with max-calls " << maxfcnUsed << " convergence for edm < " << tol + << " strategy " << options.Strategy() << std::endl; + } + + minimizer->Builder().SetPrintLevel(printLevel); + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + if (ROOT::Math::IOptions const *minuit2Opt = minuit2ExtraOptions(options)) { + int storageLevel = 1; + if (minuit2Opt->GetValue("StorageLevel", storageLevel)) { + minimizer->Builder().SetStorageLevel(storageLevel); + } + if (printLevel > 0) { + std::cout << "RooMinimizer: Minuit2 - Changing default options" << std::endl; + minuit2Opt->Print(); + } + } + + const MnStrategy strategy = makeMinuit2Strategy(options); + + _minuit2->minimum = std::make_unique(minimizer->Minimize(*_fcn, state, strategy, maxfcn, tol)); + FunctionMinimum const &minimum = *_minuit2->minimum; + + // copy minimum state (parameter values and errors) + state = minimum.UserState(); + + // Determine the status code, with the same convention as Minuit2Minimizer. + int &status = _minuit2->status; + status = 0; + std::string txt; + if (!minimum.HasPosDefCovar()) { + // this happens normally when Hesse failed + // it can happen in case MnSeed failed (see ROOT-9522) + txt = "Covar is not pos def"; + status = 5; + } + if (minimum.HasMadePosDefCovar()) { + txt = "Covar was made pos def"; + status = 1; + } + if (minimum.HesseFailed()) { + txt = "Hesse is not valid"; + status = 2; + } + if (minimum.IsAboveMaxEdm()) { + txt = "Edm is above max"; + status = 3; + } + if (minimum.HasReachedCallLimit()) { + txt = "Reached call limit"; + status = 4; + } + + MnPrint print("RooMinimizer::minimize", printLevel); + const bool validMinimum = minimum.IsValid(); + if (validMinimum) { + // print a warning message in case something is not ok + if (status != 0 && printLevel > 0) + print.Warn(txt); + } else { + // minimum is not valid when state is not valid and edm is over max or has passed call limits + if (status == 0) { + // this should not happen + txt = "unknown failure"; + status = 6; + } + print.Warn("Minimization did NOT converge,", txt); + } + + if (printLevel >= 1) { + std::cout << "RooMinimizer: Minuit2 " << (validMinimum ? "valid" : "invalid") << " minimum - status = " << status + << std::endl; + const int prec = std::cout.precision(18); + std::cout << "FVAL = " << state.Fval() << std::endl; + std::cout << "Edm = " << state.Edm() << std::endl; + std::cout.precision(prec); + std::cout << "Nfcn = " << state.NFcn() << std::endl; + if (validMinimum) { + for (MinuitParameter const &par : state.MinuitParameters()) { + std::cout << par.Name() << "\t = " << par.Value() << "\t "; + if (par.IsFixed()) + std::cout << "(fixed)" << std::endl; + else if (par.IsConst()) + std::cout << "(const)" << std::endl; + else if (par.HasLimits()) + std::cout << "+/- " << par.Error() << "\t(limited)" << std::endl; + else + std::cout << "+/- " << par.Error() << std::endl; + } + } + } + + return validMinimum; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run Minuit2's HESSE. If a function minimum from a previous minimization +/// exists, it is updated with the result. Returns whether a valid covariance +/// matrix was obtained. + +bool RooMinimizer::minuit2Hesse() +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + + if (!_minuit2->initialized) { + initMinuit2(); + } + + const unsigned int maxfcn = options.MaxFunctionCalls(); + const int printLevel = options.PrintLevel(); + + MnPrint print("RooMinimizer::hesse", printLevel); + print.Info("Using max-calls", maxfcn); + + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + MnUserParameterState &state = _minuit2->state; + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + setMinuit2ErrorDef(options.ErrorDef()); + + MnHesse hesse(makeMinuit2Strategy(options)); + + if (_minuit2->minimum) { + // run hesse and function minimum will be updated with Hesse result + hesse(*_fcn, *_minuit2->minimum, maxfcn); + state = _minuit2->minimum->UserState(); + } else { + // run Hesse on point stored in current state (independent of function minimum validity) + state = hesse(*_fcn, state, maxfcn); + } + + if (printLevel >= 3) { + std::cout << "RooMinimizer::hesse - State returned from Hesse " << std::endl; + std::cout << state << std::endl; + } + + const int covStatus = state.CovarianceStatus(); + std::string covStatusType = "not valid"; + if (covStatus == 1) + covStatusType = "approximate"; + if (covStatus == 2) + covStatusType = "full but made positive defined"; + if (covStatus == 3) + covStatusType = "accurate"; + if (covStatus == 0) + covStatusType = "full but not positive defined"; + + if (!state.HasCovariance()) { + // if false means error is not valid and this is due to a failure in Hesse + // update minimizer error status + int hstatus = 4; + // information on error state can be retrieved only if the minimum is available + if (_minuit2->minimum) { + if (_minuit2->minimum->Error().HesseFailed()) + hstatus = 1; + if (_minuit2->minimum->Error().InvertFailed()) + hstatus = 2; + else if (!(_minuit2->minimum->Error().IsPosDef())) + hstatus = 3; + } + + print.Warn("Hesse failed - matrix is", covStatusType); + print.Warn(hstatus); + + _minuit2->status += 100 * hstatus; + return false; + } + + print.Info("Hesse is valid - matrix is", covStatusType); + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run Minuit2's MINOS for the parameter with the given index. If a new +/// minimum is found, the minimization is repeated from there and MINOS is run +/// again. Returns whether both errors are valid. + +bool RooMinimizer::minuit2Minos(unsigned int index, double &errLow, double &errUp) +{ + using namespace ROOT::Minuit2; + + errLow = 0; + errUp = 0; + + const int printLevel = _config.MinimizerOptions().PrintLevel(); + MnPrint print("RooMinimizer::minos", printLevel); + + if (!_minuit2->minimum) { + print.Error("Failed - no function minimum existing"); + return false; + } + + MnUserParameterState &state = _minuit2->state; + + // need to know if parameter is const or fixed + if (isFixedOrConst(state.Parameter(index))) { + return false; + } + + if (!_minuit2->minimum->IsValid()) { + print.Error("Failed - invalid function minimum"); + return false; + } + + setMinuit2ErrorDef(_config.MinimizerOptions().ErrorDef()); + + int mstatus = runMinuit2Minos(index, errLow, errUp); + + // run again the Minimization in case of a new minimum + // bit 8 is set + if ((mstatus & 8) != 0) { + print.Info([&](std::ostream &os) { + os << "Found a new minimum: run again the Minimization starting from the new point"; + os << "\nFVAL = " << state.Fval(); + for (MinuitParameter const &par : state.MinuitParameters()) { + os << '\n' << par.Name() << "\t = " << par.Value(); + } + }); + // release parameter that was fixed in the returned state from Minos + state.Release(index); + if (!minuit2Minimize()) + return false; + // run again Minos from new Minimum (also lower error needs to be re-computed) + print.Info("Run now again Minos from the new found Minimum"); + mstatus = runMinuit2Minos(index, errLow, errUp); + + // do not reset new minimum bit to flag for other parameters + mstatus |= 8; + } + + _minuit2->status += 10 * mstatus; + _minuit2->minosStatus = mstatus; + + return ((mstatus & 1) == 0) && ((mstatus & 2) == 0); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run MINOS once for the parameter with the given index and return the MINOS +/// status bits, with the same convention as Minuit2Minimizer: +/// bit 1: lower error invalid, bit 2: upper error invalid, bit 4: invalid +/// because the maximum number of function calls was reached, bit 8: invalid +/// because a new minimum was found, bit 16: parameter is at a limit. + +int RooMinimizer::runMinuit2Minos(unsigned int index, double &errLow, double &errUp) +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + const int printLevel = options.PrintLevel(); + MnPrint print("RooMinimizer::minos", printLevel); + + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + MnUserParameterState &state = _minuit2->state; + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + // Like Minuit2Minimizer, use the default strategy for Minos. + MnMinos minos(*_fcn, *_minuit2->minimum); + + const unsigned int maxfcn = options.MaxFunctionCalls(); + // Tolerance for the migrad calls inside Minos. Cut off too small values, + // which are not needed. + const double tol = std::max(options.Tolerance(), 0.01); + + const char *parName = state.Name(index); + + if (printLevel >= 1) { + // get the real number of maxfcn used (defined in MnMinos) to be printed + unsigned int maxfcnUsed = maxfcn; + if (maxfcnUsed == 0) { + const unsigned int nvar = state.VariableParameters(); + maxfcnUsed = 2 * (nvar + 1) * (200 + 100 * nvar + 5 * nvar * nvar); + } + std::cout << "RooMinimizer::minos - Run MINOS for parameter #" << index << " : " << parName << " using max-calls " + << maxfcnUsed << ", tolerance " << tol << std::endl; + } + + const MinosError me = minos.Minos(index, maxfcn, tol); + + // Note that the only invalid condition can happen when the (npar-1) minimization fails + // The error is also invalid when the maximum number of calls is reached or a new function minimum is found + // in case of the parameter at the limit the error is not invalid. + // When the error is invalid the returned error is the Hessian error. + if (!me.LowerValid()) { + print.Warn("Invalid lower error for parameter", parName); + } + if (!me.UpperValid()) { + print.Warn("Invalid upper error for parameter", parName); + } + if (me.AtLowerLimit()) { + print.Warn("Lower error for parameter", parName, "is at the Lower limit!"); + } + if (me.AtUpperLimit()) { + print.Warn("Upper error for parameter", parName, "is at the Upper limit!"); + } + if (me.AtLowerMaxFcn()) { + print.Warn("Maximum number of function calls exceeded when running for lower error for parameter", parName); + } + if (me.AtUpperMaxFcn()) { + print.Warn("Maximum number of function calls exceeded when running for upper error for parameter", parName); + } + if (me.LowerNewMin()) { + print.Warn("New Minimum found while running Minos for lower error for parameter", parName); + } + if (me.UpperNewMin()) { + print.Warn("New Minimum found while running Minos for upper error for parameter", parName); + } + if (printLevel >= 1) { + if (me.LowerValid()) + std::cout << "Minos: Lower error for parameter " << parName << " : " << me.Lower() << std::endl; + if (me.UpperValid()) + std::cout << "Minos: Upper error for parameter " << parName << " : " << me.Upper() << std::endl; + } + + int mstatus = 0; + if (!me.LowerValid()) { + mstatus |= 1; + if (me.AtLowerMaxFcn()) + mstatus |= 4; + if (me.LowerNewMin()) + mstatus |= 8; + } + if (!me.UpperValid()) { + mstatus |= 2; + if (me.AtUpperMaxFcn()) + mstatus |= 4; + if (me.UpperNewMin()) + mstatus |= 8; + } + if (me.AtUpperLimit() || me.AtLowerLimit()) + mstatus |= 16; + + errLow = me.Lower(); + errUp = me.Upper(); + + // in case of new minimum found update also the minimum state + if (me.LowerNewMin() && me.UpperNewMin()) { + // take state with lower function value + state = (me.LowerState().Fval() < me.UpperState().Fval()) ? me.LowerState() : me.UpperState(); + } else if (me.LowerNewMin()) { + state = me.LowerState(); + } else if (me.UpperNewMin()) { + state = me.UpperState(); + } + + return mstatus; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Compute a contour for parameters ipar and jpar at the given error level +/// with Minuit2's MnContours, which requires a valid function minimum. + +bool RooMinimizer::minuit2Contour(unsigned int ipar, unsigned int jpar, unsigned int npoints, double errorDef, + double *x, double *y) +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + const int printLevel = options.PrintLevel(); + MnPrint print("RooMinimizer::contour", printLevel); + + if (!_minuit2->minimum) { + print.Error("No function minimum existing; must minimize function before"); + return false; + } + + if (!_minuit2->minimum->IsValid()) { + print.Error("Invalid function minimum"); + return false; + } + + setMinuit2ErrorDef(errorDef); + + print.Info("Computing contours at level -", errorDef); + + // switch off Minuit2 printing (for level of 0,1) + Minuit2PrintLevelRAII printLevelRAII{printLevel - 1}; + + if (options.Precision() > 0) { + _minuit2->state.SetPrecision(options.Precision()); + } + + // eventually one should specify tolerance in contours + MnContours contour(*_fcn, *_minuit2->minimum, makeMinuit2Strategy(options)); + + std::vector> result = contour(ipar, jpar, npoints); + if (result.size() != npoints) { + print.Error("Invalid result from MnContours"); + return false; + } + for (unsigned int i = 0; i < npoints; ++i) { + x[i] = result[i].first; + y[i] = result[i].second; + } + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fill the internal fit result from the Minuit2 parameter state. + +void RooMinimizer::fillResultFromMinuit2(bool isValid) +{ + using namespace ROOT::Minuit2; + + MnUserParameterState const &state = _minuit2->state; + std::vector const ¶ms = state.MinuitParameters(); + const unsigned int npar = params.size(); + + _result->fValid = isValid; + _result->fStatus = _minuit2->status; + _result->fCovStatus = minuit2CovMatrixStatus(_minuit2->minimum.get(), state); + _result->fVal = state.Fval(); + _result->fEdm = state.Edm(); + _result->fMinimType = _config.MinimizerName(); + + _result->fParams.resize(npar); + _result->fErrors.resize(npar); + for (unsigned int i = 0; i < npar; ++i) { + _result->fParams[i] = params[i].Value(); + _result->fErrors[i] = isFixedOrConst(params[i]) ? 0. : params[i].Error(); + } + + // check for fixed parameters + for (unsigned int ipar = 0; ipar < npar; ++ipar) { + if (_config.ParSettings(ipar).IsFixed()) + _result->fFixedParams[ipar] = true; + } + + // fill error matrix in packed lower-triangular storage + _result->fCovMatrix.clear(); + if (_result->fCovStatus != 0) { + _result->fCovMatrix.reserve(npar * (npar + 1) / 2); + for (unsigned int i = 0; i < npar; ++i) { + for (unsigned int j = 0; j <= i; ++j) { + _result->fCovMatrix.push_back(minuit2Covariance(state, i, j)); + } + } + } + + // global correlation coefficients, zero for fixed parameters + _result->fGlobalCC.clear(); + MnGlobalCorrelationCoeff const globalCC = state.GlobalCC(); + if (globalCC.IsValid()) { + _result->fGlobalCC.resize(npar); + for (unsigned int i = 0; i < npar; ++i) { + _result->fGlobalCC[i] = isFixedOrConst(params[i]) ? 0. : globalCC.GlobalCC()[state.IntOfExt(i)]; + } + } + // minos errors are set separately when calling calculateMinosErrors() +} + double RooMinimizer::FitResult::lowerError(unsigned int i) const { // return lower Minos error for parameter i diff --git a/roofit/roofitcore/src/RooMinimizerFcn.cxx b/roofit/roofitcore/src/RooMinimizerFcn.cxx index 6a3235950f2b3..3f482f65b515f 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooMinimizerFcn.cxx @@ -16,8 +16,10 @@ ////////////////////////////////////////////////////////////////////////////// /// \class RooMinimizerFcn -/// RooMinimizerFcn is an interface to the ROOT::Math::IBaseFunctionMultiDim, -/// a function that ROOT's minimisers use to carry out minimisations. +/// RooMinimizerFcn adapts a RooAbsReal to the function interface of the +/// minimizers used by RooMinimizer: the ROOT::Minuit2::FCNBase interface for +/// Minuit2, and ROOT::Math::IBaseFunctionMultiDim for the other minimizers +/// that are reached via ROOT::Math::Minimizer. /// #include "RooMinimizerFcn.h" @@ -35,7 +37,6 @@ #include "RooRealVar.h" #include "Math/Functor.h" -#include "Minuit2/Minuit2Minimizer.h" #include "TMatrixDSym.h" #include @@ -151,14 +152,18 @@ RooMinimizerFcn::RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context) { unsigned int nDim = getNDim(); - if (context->_cfg.useGradient && funct->hasGradient()) { + _useGradient = context->_cfg.useGradient && funct->hasGradient(); + _useHessian = context->_cfg.useHessian && funct->hasHessian(); + + auto evalFunc = [this](const double *x) { return (*this)(x); }; + if (_useGradient) { _gradientOutput.resize(_allParams.size()); - _multiGenFcn = std::make_unique(this, &RooMinimizerFcn::operator(), - &RooMinimizerFcn::evaluateGradient, nDim); + auto gradFunc = [this](const double *x, double *out) { evaluateGradient(x, out); }; + _multiGenFcn = std::make_unique(evalFunc, nDim, gradFunc); } else { - _multiGenFcn = std::make_unique(std::cref(*this), nDim); + _multiGenFcn = std::make_unique(evalFunc, nDim); } - if (context->_cfg.useHessian) { + if (_useHessian) { _hessianOutput.resize(_allParams.size() * _allParams.size()); } } @@ -267,7 +272,7 @@ RooArgSet RooMinimizerFcn::freezeDisconnectedParameters() const return changedSet; } -bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) const +std::vector RooMinimizerFcn::Hessian(std::vector const &x) const { // Set the parameter values for this iteration for (unsigned index = 0; index < getNDim(); index++) { @@ -280,6 +285,7 @@ bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) co std::size_t m = _allParamsInit.size(); std::size_t n = getNDim(); + std::vector out(n * n); std::size_t iAll = 0; std::size_t iFloating = 0; for (RooAbsArg *param_i : _allParamsInit) { @@ -308,24 +314,19 @@ bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) co std::cout << std::endl; } } - return true; + return out; +} + +std::vector RooMinimizerFcn::Gradient(std::vector const &x) const +{ + std::vector out(getNDim()); + evaluateGradient(x.data(), out.data()); + return out; } -void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer *context) +void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim) const { minim.SetFunction(*_multiGenFcn); - if (context->_cfg.useHessian && _funct->hasHessian()) { - minim.SetHessianFunction( - std::bind(&RooMinimizerFcn::evaluateHessian, this, std::placeholders::_1, std::placeholders::_2)); - } - // The independence information for skipping vanishing second derivatives - // in numerical Hessian computations is a Minuit2-only feature, so it is - // wired up directly with the concrete minimizer type instead of going - // through the ROOT::Math::Minimizer interface. - if (auto *minuit2 = dynamic_cast(&minim)) { - minuit2->SetSecondDerivativeAlwaysVanishesFunc( - [this](unsigned int i, unsigned int j) { return secondDerivativeAlwaysVanishes(i, j); }); - } } //////////////////////////////////////////////////////////////////////////////// @@ -365,7 +366,7 @@ void RooMinimizerFcn::buildSecondDerivMask() const /// (indices in the space of all floatable parameters, matching Minuit's /// external parameter indices) is identically zero because the parameters /// share no additive term of the minimized function. -bool RooMinimizerFcn::secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const +bool RooMinimizerFcn::SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const { std::call_once(_secondDerivMaskOnce, &RooMinimizerFcn::buildSecondDerivMask, this); return !_secondDerivMask[getNDim() * i + j]; diff --git a/roofit/roofitcore/src/RooMinimizerFcn.h b/roofit/roofitcore/src/RooMinimizerFcn.h index de7052c161198..760a9015d9abb 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.h +++ b/roofit/roofitcore/src/RooMinimizerFcn.h @@ -22,6 +22,8 @@ #include "RooAbsReal.h" #include "RooArgList.h" +#include + #include #include #include @@ -36,7 +38,10 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { public: RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context); - void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) override; + /// Set this function on a ROOT::Math::Minimizer. Only used for minimizers + /// other than Minuit2, which RooMinimizer drives directly via the + /// ROOT::Minuit2::FCNBase interface. + void initMinimizer(ROOT::Math::Minimizer &minim) const; std::string getFunctionName() const override; std::string getFunctionTitle() const override; @@ -45,16 +50,26 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { double operator()(const double *x) const; void evaluateGradient(const double *x, double *out) const; - bool evaluateHessian(std::span x, double *out) const; RooArgSet freezeDisconnectedParameters() const override; - bool secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double operator()(std::vector const &x) const override { return (*this)(x.data()); } + bool HasGradient() const override { return _useGradient; } + std::vector Gradient(std::vector const &x) const override; + bool HasHessian() const override { return _useHessian; } + std::vector Hessian(std::vector const &x) const override; + bool SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const override; + /// @} private: void buildSecondDerivMask() const; RooAbsReal *_funct = nullptr; + bool _useGradient = false; ///< Whether to provide the analytic gradient to the minimizer. + bool _useHessian = false; ///< Whether to provide the analytic Hessian to the minimizer. + /// Adapter to ROOT::Math::Minimizer, only used for minimizers other than Minuit2. std::unique_ptr _multiGenFcn; mutable std::vector _gradientOutput; mutable std::vector _hessianOutput; diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index 37b69e3acf2f3..242d37e72482d 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -21,7 +21,7 @@ #include "RooMsgService.h" #include "RooMinimizer.h" -#include "Minuit2/Minuit2Minimizer.h" +#include "../RooAbsMinimizerFcn.h" #include "Minuit2/MnStrategy.h" #include @@ -169,9 +169,8 @@ void LikelihoodGradientJob::update_state() auto gradient_message = get_manager()->messenger().receive_from_master_on_worker(&more); assert(more); - auto gradient_message_begin = gradient_message.data(); - auto gradient_message_end = - gradient_message_begin + gradient_message.size() / sizeof(ROOT::Minuit2::DerivatorElement); + auto gradient_message_begin = gradient_message.data(); + auto gradient_message_end = gradient_message_begin + gradient_message.size() / sizeof(DerivatorElement); std::copy(gradient_message_begin, gradient_message_end, grad_.begin()); auto minuit_internal_x_message = @@ -193,9 +192,6 @@ void LikelihoodGradientJob::update_state() std::copy(offsets_message_begin, offsets_message_end, shared_offset_.offsets().begin()); } - // Since the gradient parallelization only support Minuit 2, we can do this cast - auto &minim = static_cast(*minimizer_->_minimizer); - // The master already knows the function value at the current point from // the line search that Minuit just completed; pre-seeding the derivator's // cache with it makes the SetupDifferentiate call below skip its full @@ -206,8 +202,8 @@ void LikelihoodGradientJob::update_state() // note: the next call must stay after the (possible) update of the offset, because it // calls the likelihood function, so the offset must be correct at this point - gradf_.SetupDifferentiate(minimizer_->getNPar(), minim.GetFCN(), minuit_internal_x_.data(), - minimizer_->fitter()->Config().ParamsSettings()); + gradf_.SetupDifferentiate(minimizer_->getNPar(), minimizer_->_fcn.get(), minuit_internal_x_.data(), + minimizer_->_config.ParamsSettings()); } } @@ -218,11 +214,8 @@ void LikelihoodGradientJob::update_state() void LikelihoodGradientJob::run_derivator(unsigned int i_component) const { - // Since the gradient parallelization only support Minuit 2, we can do this cast - auto &minim = static_cast(*minimizer_->_minimizer); - // Calculate the derivative etc for these parameters - grad_[i_component] = gradf_.FastPartialDerivative(minim.GetFCN(), minimizer_->fitter()->Config().ParamsSettings(), + grad_[i_component] = gradf_.FastPartialDerivative(minimizer_->_fcn.get(), minimizer_->_config.ParamsSettings(), i_component, grad_[i_component]); } diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h index da94b01efb800..5951d64b8eb8d 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h @@ -17,9 +17,10 @@ #include "RooFit/TestStatistics/LikelihoodGradientWrapper.h" #include "Math/MinimizerOptions.h" -#include "Minuit2/NumericalDerivator.h" #include "Minuit2/MnMatrix.h" +#include "NumericalDerivator.h" + #include #include @@ -63,7 +64,7 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien struct task_result_t { std::size_t job_id; std::size_t task_id; - ROOT::Minuit2::DerivatorElement grad; + DerivatorElement grad; }; void send_back_task_result_from_worker(std::size_t task) override; bool receive_task_result_on_master(const RooFit::MultiProcess::Message &message) override; @@ -75,8 +76,8 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien // members // mutables below are because ROOT::Math::IMultiGradFunction::DoDerivative is const - mutable std::vector grad_; - mutable ROOT::Minuit2::NumericalDerivator gradf_; + mutable std::vector grad_; + mutable NumericalDerivator gradf_; std::size_t N_tasks_ = 0; std::size_t N_tasks_at_workers_ = 0; diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx index d785bfe4eabce..1fdfc890da782 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -17,57 +17,11 @@ #include "RooAbsPdf.h" #include "RooNaNPacker.h" -#include -#include - #include // std::setprecision namespace RooFit { namespace TestStatistics { -namespace { - -class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { -public: - MinuitGradFunctor(MinuitFcnGrad const &fcn, double errorLevel) : _fcn{fcn}, _up{errorLevel} {} - - double operator()(std::vector const &v) const override { return _fcn(v.data()); } - double Up() const override { return _up; } - void SetErrorDef(double val) override { _up = val; } - bool HasGradient() const override { return true; } - std::vector Gradient(std::vector const ¶ms) const override - { - std::vector grad(_fcn.getNDim()); - _fcn.Gradient(params.data(), grad.data()); - return grad; - } - // Unhide the 4-argument overload from FCNBase, which forwards to Gradient(). - // Otherwise GCC's -Woverloaded-virtual (enabled with -Werror on some CI - // targets) complains because we only override the 5-argument overload here. - using FCNBase::GradientWithPrevResult; - std::vector GradientWithPrevResult(std::vector const &v, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtV) const override - { - std::vector output(v.size()); - _fcn.GradientWithPrevResult(v.data(), output.data(), previous_grad, previous_g2, previous_gstep, fValAtV); - return output; - } - ROOT::Minuit2::GradientParameterSpace gradParameterSpace() const override - { - return _fcn.returnsInMinuit2ParameterSpace() ? ROOT::Minuit2::GradientParameterSpace::Internal - : ROOT::Minuit2::GradientParameterSpace::External; - } - - // TODO: Implement this - bool SecondDerivativeAlwaysVanishes(unsigned int /*i*/, unsigned int /*j*/) const override { return false; } - -private: - MinuitFcnGrad const &_fcn; - double _up; -}; - -} // namespace - /** \class MinuitFcnGrad * * \brief Minuit-RooMinimizer interface which synchronizes parameter data and coordinates evaluation of likelihood @@ -77,8 +31,8 @@ class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { * calculator classes to evaluate likelihood and likelihood gradient values and returns them to Minuit. The Wrapper * objects do the actual calculations. These are constructed inside the MinuitFcnGrad constructor using the RooAbsL * likelihood passed in to the constructor, usually directly from RooMinimizer, with which this class is intimately - * coupled, being a RooAbsMinimizerFcn implementation. MinuitFcnGrad inherits from ROOT::Math::IMultiGradFunction as - * well, which allows it to be used as the FCN and GRAD parameters Minuit expects. + * coupled, being a RooAbsMinimizerFcn implementation. Like all RooAbsMinimizerFcn classes, MinuitFcnGrad implements + * the ROOT::Minuit2::FCNBase interface, which allows it to be used as the FCN and GRAD parameters Minuit expects. * * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). @@ -141,9 +95,9 @@ void MinuitFcnGrad::syncOffsets() const } } -double MinuitFcnGrad::operator()(const double *x) const +double MinuitFcnGrad::operator()(std::vector const &x) const { - syncParameterValuesFromMinuitCalls(x, false); + syncParameterValuesFromMinuitCalls(x.data(), false); syncOffsets(); @@ -253,23 +207,28 @@ bool MinuitFcnGrad::syncParameterValuesFromMinuitCalls(const double *x, bool min return aParamWasUpdated; } -void MinuitFcnGrad::Gradient(const double *x, double *grad) const +std::vector MinuitFcnGrad::Gradient(std::vector const &x) const { + std::vector grad(getNDim()); _calculatingGradient = true; - syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); + syncParameterValuesFromMinuitCalls(x.data(), returnsInMinuit2ParameterSpace()); syncOffsets(); - _gradient->fillGradient(grad); + _gradient->fillGradient(grad.data()); _calculatingGradient = false; + return grad; } -void MinuitFcnGrad::GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtX) const +std::vector MinuitFcnGrad::GradientWithPrevResult(std::vector const &x, double *previous_grad, + double *previous_g2, double *previous_gstep, + double fValAtX) const { + std::vector grad(getNDim()); _calculatingGradient = true; - syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); + syncParameterValuesFromMinuitCalls(x.data(), returnsInMinuit2ParameterSpace()); syncOffsets(); - _gradient->fillGradientWithPrevResult(grad, previous_grad, previous_g2, previous_gstep, fValAtX); + _gradient->fillGradientWithPrevResult(grad.data(), previous_grad, previous_g2, previous_gstep, fValAtX); _calculatingGradient = false; + return grad; } bool MinuitFcnGrad::Synchronize(std::vector ¶meters) @@ -278,15 +237,15 @@ bool MinuitFcnGrad::Synchronize(std::vector ¶m applyToLikelihood([&](auto &l) { l.synchronizeParameterSettings(parameters); }); _gradient->synchronizeParameterSettings(parameters); - applyToLikelihood([&](auto &l) { l.synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); }); - _gradient->synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); + applyToLikelihood([&](auto &l) { l.synchronizeWithMinimizer(minimizerOptions()); }); + _gradient->synchronizeWithMinimizer(minimizerOptions()); return returnee; } -void MinuitFcnGrad::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer * /*context*/) +ROOT::Minuit2::GradientParameterSpace MinuitFcnGrad::gradParameterSpace() const { - auto &minuit = dynamic_cast(minim); - minuit.SetFCN(getNDim(), std::make_unique(*this, minim.ErrorDef())); + return returnsInMinuit2ParameterSpace() ? ROOT::Minuit2::GradientParameterSpace::Internal + : ROOT::Minuit2::GradientParameterSpace::External; } } // namespace TestStatistics diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h index f345fdf67e188..7e62676456ce8 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h @@ -21,7 +21,6 @@ #include "../RooAbsMinimizerFcn.h" #include -#include "Math/IFunction.h" // ROOT::Math::IMultiGradFunction class RooMinimizer; @@ -34,19 +33,24 @@ class MinuitFcnGrad : public RooAbsMinimizerFcn { std::vector ¶meters, LikelihoodMode likelihoodMode, LikelihoodGradientMode likelihoodGradientMode); - void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) override; - /// Overridden from RooAbsMinimizerFcn to include gradient strategy synchronization. bool Synchronize(std::vector ¶meter_settings) override; bool returnsInMinuit2ParameterSpace() const { return _gradient->usesMinuitInternalValues(); } - double operator()(const double *x) const; - - /// IMultiGradFunction overrides necessary for Minuit - void Gradient(const double *x, double *grad) const; - void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtX) const; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double operator()(std::vector const &x) const override; + bool HasGradient() const override { return true; } + std::vector Gradient(std::vector const &x) const override; + // Unhide the 4-argument overload from FCNBase, which forwards to Gradient(). + // Otherwise GCC's -Woverloaded-virtual (enabled with -Werror on some CI + // targets) complains because we only override the 5-argument overload here. + using FCNBase::GradientWithPrevResult; + std::vector GradientWithPrevResult(std::vector const &x, double *previous_grad, double *previous_g2, + double *previous_gstep, double fValAtX) const override; + ROOT::Minuit2::GradientParameterSpace gradParameterSpace() const override; + /// @} inline std::string getFunctionName() const override { return _likelihood->GetName(); } diff --git a/math/minuit2/src/NumericalDerivator.cxx b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx similarity index 97% rename from math/minuit2/src/NumericalDerivator.cxx rename to roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx index d028ed95d7c84..5151cf44116bd 100644 --- a/math/minuit2/src/NumericalDerivator.cxx +++ b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx @@ -1,4 +1,3 @@ -// @(#)root/mathcore:$Id$ // Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 // /********************************************************************** @@ -23,18 +22,20 @@ * Modified things (w.r.t. original) are indicated by MODIFIED. */ -#include "Minuit2/NumericalDerivator.h" -#include -#include -#include +#include "NumericalDerivator.h" + +#include #include + +#include #include -#include "Fit/ParameterSettings.h" +#include +#include -#include // needed here because in Fitter is only a forward declaration +namespace RooFit { +namespace TestStatistics { -namespace ROOT { -namespace Minuit2 { +using ROOT::Minuit2::FCNBase; NumericalDerivator::NumericalDerivator(bool always_exactly_mimic_minuit2) : fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) @@ -286,5 +287,5 @@ std::ostream &operator<<(std::ostream &out, const DerivatorElement &value) << ", step_size: " << value.step_size << ")"; } -} // namespace Minuit2 -} // namespace ROOT +} // namespace TestStatistics +} // namespace RooFit diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.h similarity index 83% rename from math/minuit2/inc/Minuit2/NumericalDerivator.h rename to roofit/roofitcore/src/TestStatistics/NumericalDerivator.h index f8e42fc35fd7a..ad2c83001781a 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.h @@ -1,4 +1,3 @@ -// @(#)root/mathcore:$Id$ // Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 // /********************************************************************** @@ -13,10 +12,13 @@ * Authors: L. Moneta, J. T. Offermann * Modified version created on: Sep 27, 2017 * Author: E. G. P. Bos + * + * Lived in Minuit2 as ROOT::Minuit2::NumericalDerivator until it was moved + * to RooFit, its only user. */ -#ifndef ROOT_Minuit2_NumericalDerivator -#define ROOT_Minuit2_NumericalDerivator +#ifndef ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator +#define ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator #include "Fit/ParameterSettings.h" #include "Minuit2/MnParameterTransformation.h" @@ -27,8 +29,8 @@ #include -namespace ROOT { -namespace Minuit2 { +namespace RooFit { +namespace TestStatistics { // Holds all necessary derivatives and associated numbers (per parameter) used in the NumericalDerivator class. struct DerivatorElement { @@ -44,7 +46,7 @@ class NumericalDerivator { NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, bool always_exactly_mimic_minuit2 = true); - void SetupDifferentiate(unsigned int nDim, const FCNBase *function, const double *cx, + void SetupDifferentiate(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *cx, std::span parameters); /// Pre-seed the cache of the function value at the central point, so that a @@ -58,17 +60,18 @@ class NumericalDerivator { fVal = fval; } - std::vector Differentiate(unsigned int nDim, const FCNBase *function, const double *x, + std::vector Differentiate(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, + const double *x, std::span parameters, std::span previous_gradient); - DerivatorElement PartialDerivative(unsigned int nDim, const FCNBase *function, const double *x, + DerivatorElement PartialDerivative(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *x, std::span parameters, unsigned int i_component, DerivatorElement previous); - DerivatorElement FastPartialDerivative(const FCNBase *function, + DerivatorElement FastPartialDerivative(const ROOT::Minuit2::FCNBase *function, std::span parameters, unsigned int i_component, const DerivatorElement &previous); - DerivatorElement operator()(unsigned int nDim, const FCNBase *function, const double *x, + DerivatorElement operator()(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *x, std::span parameters, unsigned int i_component, const DerivatorElement &previous); @@ -114,7 +117,7 @@ class NumericalDerivator { std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); -} // namespace Minuit2 -} // namespace ROOT +} // namespace TestStatistics +} // namespace RooFit -#endif // ROOT_Minuit2_NumericalDerivator +#endif // ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator diff --git a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx index db0641eba7977..5f1f44785ae3f 100644 --- a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx +++ b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx @@ -968,10 +968,10 @@ TEST(MinuitFcnGrad, DISABLED_CompareToRooMinimizerFcn) RFTS::LikelihoodGradientMode::multiprocess); RooMinimizerFcn vanilla_fcn(nll_vanilla.get(), &m_vanilla); - EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn).data()), modularL_fcn(getParamVals(modularL_fcn).data())); + EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn)), modularL_fcn(getParamVals(modularL_fcn))); // let's also check with absolutely certain same parameter values, both of them - EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn).data()), modularL_fcn(getParamVals(vanilla_fcn).data())); - EXPECT_EQ(vanilla_fcn(getParamVals(modularL_fcn).data()), modularL_fcn(getParamVals(modularL_fcn).data())); + EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn)), modularL_fcn(getParamVals(vanilla_fcn))); + EXPECT_EQ(vanilla_fcn(getParamVals(modularL_fcn)), modularL_fcn(getParamVals(modularL_fcn))); // reset static variables to automatic RFMP::Config::LikelihoodJob::defaultNEventTasks = RFMP::Config::LikelihoodJob::automaticNEventTasks; diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 97d1b8d416199..5decb10285dc1 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -78,7 +78,9 @@ void randomizeParameters(const RooArgSet ¶meters) double mul = rng.Uniform(lowerBound, upperBound); auto par = dynamic_cast(param); - if (!par) + // Constant parameters must not be varied. Note that they can also have + // an infinite range, so the randomization logic below would break. + if (!par || par->isConstant()) continue; double val = par->getVal(); val = val + mul * (mul > 0 ? (par->getMax() - val) : (val - par->getMin())); @@ -97,12 +99,15 @@ class FactoryTestParams { public: FactoryTestParams() = default; FactoryTestParams(std::string const &name, WorkspaceSetupFunc setupWorkspace, CreateNLLFunc createNLL, - double fitResultTolerance, bool randomizeParameters, std::size_t nEvents = 100) + double fitResultTolerance, bool randomizeParameters, double hesseTolerance = -1., + bool enableHessian = true, std::size_t nEvents = 100) : _name{name}, _setupWorkspace{setupWorkspace}, _createNLL{createNLL}, _fitResultTolerance{fitResultTolerance}, + _hesseTolerance{hesseTolerance}, _randomizeParameters{randomizeParameters}, + _enableHessian{enableHessian}, _nEvents{nEvents} { } @@ -111,7 +116,15 @@ class FactoryTestParams { WorkspaceSetupFunc _setupWorkspace; CreateNLLFunc _createNLL; double _fitResultTolerance = 1e-4; + // Tolerance for comparing parameter errors, which are less precise than + // the values because Minuits numeric Hesse is only accurate to about the + // percent level. If negative, _fitResultTolerance is used. + double _hesseTolerance = -1.; bool _randomizeParameters = true; + // Whether to test the minimization also with the analytic Hessian from + // Clad. Disabled for models where the Hessian is known to come out wrong + // or can't be generated at all (see comments at the test definitions). + bool _enableHessian = true; std::size_t _nEvents = 100; }; @@ -131,14 +144,16 @@ class FactoryTest : public testing::TestWithParam { std::unique_ptr _changeMsgLvl; }; -std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useGradient = true) +std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useAD = true, bool useHessian = true) { RooMinimizer::Config cfg; - cfg.useGradient = useGradient; + cfg.useGradient = useAD; + cfg.useHessian = useAD && useHessian; RooMinimizer m{absReal, cfg}; m.setPrintLevel(-1); m.setStrategy(0); m.minimize("Minuit2"); + m.hesse(); // to use the analytical Hessian for error estimation return std::unique_ptr{m.save()}; } @@ -166,8 +181,13 @@ TEST_P(FactoryTest, NLLFit) // We want to use the generated code also for the nominal likelihood. Like // this, we make sure to validate also the NLL values of the generated code. static_cast(*nllFunc).setUseGeneratedFunctionCode(true); + // Let's also validate the Hessian + if (_params._enableHessian) { + static_cast(*nllFunc).generateHessian(); + } double tol = _params._fitResultTolerance; + double hesseTol = _params._hesseTolerance < 0. ? tol : _params._hesseTolerance; EXPECT_NEAR(nllRef->getVal(observables), nllFunc->getVal(), tol); @@ -204,7 +224,7 @@ TEST_P(FactoryTest, NLLFit) paramsRefNll.assign(parametersOrig); // Minimize the RooFuncWrapper Implementation with AD - auto resultAd = runMinimizer(*nllFunc); + auto resultAd = runMinimizer(*nllFunc, true, _params._enableHessian); paramsRefNll.assign(parametersOrig); // Minimize the reference NLL @@ -216,8 +236,8 @@ TEST_P(FactoryTest, NLLFit) // because for very small correlations it's usually not the same within the // relative tolerance because you would compare two small values that are // only different from zero because of noise. - EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, tol)); - EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, tol)); + EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, hesseTol)); + EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, hesseTol)); } /// Initial minimization that was not based on any other tutorial/test. @@ -225,7 +245,10 @@ FactoryTestParams param1{"Gaussian", [](RooWorkspace &ws) { constexpr double inf = std::numeric_limits::infinity(); ws.import(RooRealVar{"mu", "mu", -inf, inf}); - ws.factory("sum::mu_shifted(mu, shift[1.0, -10, 10])"); + // The shift parameter needs to be constant, because otherwise the + // model would only depend on the sum "mu + shift": the Hessian at the + // minimum is singular and the parameter errors are meaningless. + ws.factory("sum::mu_shifted(mu, shift[1.0])"); ws.factory("prod::sigma_scaled(sigma[3.0, 0.01, 10], 1.5)"); ws.factory("Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"); @@ -236,7 +259,8 @@ FactoryTestParams param1{"Gaussian", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/false}; + /*randomizeParameters=*/false, + /*hesseTolerance=*/1e-3}; /// Test based on the rf301 tutorial. FactoryTestParams param2{"PolyVar", @@ -271,7 +295,8 @@ FactoryTestParams param4{"ConstraintSum", pdf.createNLL(data, ExternalConstraints(*ws.pdf("fconstext")), backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-3}; namespace { @@ -343,12 +368,14 @@ void getSimPdfModel(RooWorkspace &ws) } // namespace /// Test based on the simultaneous fit shown in CHEP'23 results -FactoryTestParams param5{"SimPdf", getSimPdfModel, +FactoryTestParams param5{"SimPdf", + getSimPdfModel, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 5e-3, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/5e-2}; FactoryTestParams param6{"GaussianExtended", [](RooWorkspace &ws) { @@ -390,12 +417,14 @@ void getDataHistModel(RooWorkspace &ws) } // namespace /// Test based on rf706 tutorial -FactoryTestParams param7{"HistPdf", getDataHistModel, +FactoryTestParams param7{"HistPdf", + getDataHistModel, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/2e-3}; FactoryTestParams param8{"Lognormal", [](RooWorkspace &ws) { @@ -406,7 +435,8 @@ FactoryTestParams param8{"Lognormal", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams param8p1{"LognormalStandard", [](RooWorkspace &ws) { @@ -418,7 +448,8 @@ FactoryTestParams param8p1{"LognormalStandard", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 3e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/3e-2}; FactoryTestParams param9{"Poisson", [](RooWorkspace &ws) { @@ -429,9 +460,10 @@ FactoryTestParams param9{"Poisson", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; -// A RooPoisson where x is not rounded, like it is used in HistFactory +// A RooPoisson where x is not rounded, like it is used in HistFactory. FactoryTestParams param10{"PoissonNoRounding", [](RooWorkspace &ws) { constexpr double inf = std::numeric_limits::infinity(); @@ -446,7 +478,8 @@ FactoryTestParams param10{"PoissonNoRounding", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams param11{"ClassFactory1D", [](RooWorkspace &ws) { @@ -505,10 +538,14 @@ void func_1_grad(double const * /*x*/, double *grad_out) grad_out[0] = 1; grad_out[1] = -1; } +void func_1_hess(double const * /*x*/, double *hess_out) +{ + std::fill(hess_out, hess_out + 4, 0.0); +} auto &functor_1() { - static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad}; + static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad, func_1_hess}; return functor; } @@ -532,9 +569,35 @@ void func_gaussian_grad(double const *x, double *grad_out) grad_out[2] = f * arg * arg * inv_sig2 / sig; } +void func_gaussian_hess(double const *x, double *hess_out) +{ + const double arg = x[0] - x[1]; + const double sig = x[2]; + + const double invSig2 = 1.0 / (sig * sig); + const double t = arg * arg * invSig2; + const double f = std::exp(-0.5 * t); + + const double h00 = f * (t - 1.0) * invSig2; + const double h02 = f * arg * (2.0 - t) * invSig2 / sig; + const double h22 = f * t * (t - 3.0) * invSig2; + + // The function is invariant under exchanging x[0] - x[1] with x[1] - x[0], + // so the Hessian rows/columns for x[0] and x[1] only differ in sign. + hess_out[0] = h00; + hess_out[1] = -h00; + hess_out[2] = h02; + hess_out[3] = -h00; + hess_out[4] = h00; + hess_out[5] = -h02; + hess_out[6] = h02; + hess_out[7] = -h02; + hess_out[8] = h22; +} + auto &functor_gaussian() { - static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad}; + static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad, func_gaussian_hess}; return functor; } @@ -546,10 +609,14 @@ double func_1_1D_diff(double /*x*/) { return 1.; } +double func_1_1D_diff2(double /*x*/) +{ + return 0.; +} auto &functor_1_1D() { - static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff}; + static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff, func_1_1D_diff2}; return functor; } @@ -558,12 +625,17 @@ auto &functor_1_1D() FactoryTestParams param13{"RooFunctor", [](RooWorkspace &ws) { RooRealVar x("x", "", 0.0, -4, 4); - RooRealVar mu("mu", "", 0.0, -4, 4); + // The initial mu value is away from zero such that relative + // comparisons of the fitted values are meaningful. + RooRealVar mu("mu", "", 2.0, -4, 4); RooRealVar shift("shift", "", 1.0, -4, 4); shift.setConstant(true); RooFunctorBinding mu_shifted("mu_shifted", "", functor_1(), {mu, shift}); RooFunctor1DBinding mu_shifted_1D("mu_shifted_1D", "", functor_1_1D(), {mu}); - RooRealVar sigma("sigma", "", 4., 0.01, 10.); + // The initial sigma value is small enough for the Gaussian to be + // well-contained in the observable range, such that the generated + // data constrains all parameters and the fits are well-defined. + RooRealVar sigma("sigma", "", 1.5, 0.01, 10.); RooFunctorPdfBinding gauss_1("model_1", "", functor_gaussian(), {x, mu_shifted, sigma}); RooFunctorPdfBinding gauss_2("model_2", "", functor_gaussian(), {x, mu_shifted_1D, sigma}); @@ -574,15 +646,27 @@ FactoryTestParams param13{"RooFunctor", ws.import(model); ws.defineSet("observables", vars); + + // The model is statistically very flat around the minimum, so with + // only the default 100 events, the fitted values with and without + // AD can differ at the percent level just because the minimizers + // stop at slightly different points in the flat valley. Generate + // more events to make the minimum well-defined. + std::unique_ptr data{ws.pdf("model")->generate(*ws.set("observables"), 1000)}; + std::unique_ptr binned{data->binnedClone()}; + binned->SetName("data"); + ws.import(*binned); }, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-3, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams makeTestParams(const char *name, std::vector const &expressions, - double fitResultTolerance, bool randomizeParameters = true, std::size_t nEvents = 100) + double fitResultTolerance, bool randomizeParameters = true, + double hesseTolerance = -1., bool enableHessian = true, std::size_t nEvents = 100) { return {name, [=](RooWorkspace &ws) { @@ -597,6 +681,8 @@ FactoryTestParams makeTestParams(const char *name, std::vector cons }, fitResultTolerance, randomizeParameters, + hesseTolerance, + enableHessian, nEvents}; } @@ -605,12 +691,15 @@ auto testValues = testing::Values( makeTestParams( "BifurGauss", {"x[0, -10, 10]", "mu[0, -10, 10]", "BifurGauss::model(x, mu, sigmaL[3.0, 0.01, 10], sigmaR[2.0, 0.01, 10])"}, - 1e-4, false), + 1e-4, false, /*hesseTolerance=*/2e-2), + // The shift parameter needs to be constant, because otherwise the model + // would only depend on the sum "mu + shift": the Hessian at the minimum is + // singular and the parameter errors are meaningless. makeTestParams("RooFormulaVar", {"expr::mu_shifted('mu+shift',{mu[0, -10, 10], shift[1.0]})", "expr::sigma_scaled('sigma*1.5',{sigma[3.0, 0.01, 10]})", "Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"}, - 1e-4, false), + 1e-4, false, /*hesseTolerance=*/5e-3), // Test for the uniform pdf. Since it doesn't depend on any parameters, we need // to add it to some other model like a Gaussian to get a meaningful fit model. makeTestParams("Uniform", @@ -624,28 +713,39 @@ auto testValues = testing::Values( "Gaussian::sig4(x, 6.0, sigma4[2.0, .01, 10])", "RecursiveFraction::recfrac({a1[0.25, 0.0, 1.0], a2[0.25, 0.0, 1.0]})", "SUM::model(a1 * sig1, a2 * sig2, recfrac * sig3, sig4)"}, - 5e-3, true), + 5e-3, true, /*hesseTolerance=*/3e-2), makeTestParams("RooCBShape", {"x[0., -200., 200.]", "x0[100., -200., 200.]", "CBShape::model(x, x0, sigma[2., 1.E-1, 100.], alpha[1., 1.E-1, 100.], n[1., 1.E-1, 100.])"}, 6e-3, true), - makeTestParams("RooBernstein", - {"Bernstein::model(x[0., 100.], {c0[0.3, 0., 10.], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"}, 6e-3, - true), + // The first Bernstein coefficient needs to be constant, because the pdf is + // invariant under a common rescaling of all coefficients: the Hessian at + // the minimum is singular and the parameter errors are meaningless. + makeTestParams("RooBernstein", {"Bernstein::model(x[0., 100.], {c0[0.3], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"}, + 6e-3, true, + /*hesseTolerance=*/2e-2), // We're testing several Landau configurations, because the underlying // ROOT::Math::landau_cdf is defined piecewise. Like this, we're covering // all possible code paths in the pullback. - makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false), - makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false), - // With sl = 10, the Landau is nearly uniform over the observable window, so - // with only 100 events the fit is unidentifiable: the minimizers stop at - // arbitrary points in an almost flat valley, and the fitted values can't be - // meaningfully compared. Generate more events to make the minimum - // well-defined. + makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + // With sl = 10, the Landau is nearly uniform over the observable window, + // so with only 100 events the fit is unidentifiable: the minimizers stop + // at arbitrary points in an almost flat valley, and the fitted values + // can't be meaningfully compared. Generate more events to make the + // minimum well-defined. + // The Hesse tolerance is wider than for the other Landau tests: with sl = + // 10, the model is poorly conditioned and Minuits numeric Hesse in the + // reference fit is less accurate (the Clad Hessian was validated against + // finite differences of the exact gradient). makeTestParams("RooLandau3", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[10., 0.01, 50.])"}, 7e-3, false, - 10000), - makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false), - makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false), + /*hesseTolerance=*/3e-2, /*enableHessian=*/true, /*nEvents=*/10000), + makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), makeTestParams( "RooRealSumPdf1", {"Gaussian::gx(x[-10,10],m[0],1.0)", "Chebychev::ch(x,{0.1,0.2,-0.3})", "RealSumPdf::model({gx, ch}, {f[0,1]})"}, diff --git a/roofit/roofitcore/test/testRooMinimizer.cxx b/roofit/roofitcore/test/testRooMinimizer.cxx index f3f3d3cbecee1..aeaf78895987d 100644 --- a/roofit/roofitcore/test/testRooMinimizer.cxx +++ b/roofit/roofitcore/test/testRooMinimizer.cxx @@ -159,7 +159,7 @@ INSTANTIATE_TEST_SUITE_P(RooMinimizer, EvalBackendParametrizedTest, testing::Val // parameter is fixed after the RooMinimizer was constructed. The latter case // is a regression test for the translation between Minuit-internal parameter // indices (which exclude fixed parameters) and the external indices that -// RooMinimizerFcn::secondDerivativeAlwaysVanishes() is defined in. +// RooMinimizerFcn::SecondDerivativeAlwaysVanishes() is defined in. TEST(RooMinimizer, SecondDerivativeAlwaysVanishesHesse) { RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; diff --git a/roofit/xroofit/src/xRooFit.cxx b/roofit/xroofit/src/xRooFit.cxx index 6a526d549b289..f6da742d6e298 100644 --- a/roofit/xroofit/src/xRooFit.cxx +++ b/roofit/xroofit/src/xRooFit.cxx @@ -1185,7 +1185,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, }*/ // std::cout << "nIterations = " << _minimizer->fitter()->GetMinimizer()->NIterations() << std::endl; - // std::cout << "covQual before hesse = " << _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() << + // std::cout << "covQual before hesse = " << _minimizer->fitter()->Result().CovMatrixStatus() << // std::endl; sIdx = -1; if (hesseStrategy == -1) { @@ -1204,7 +1204,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (strategy == 2 && hesseStrategy == 2) { // don't repeat hesse if strategy=2 and hesseStrategy=2, and the matrix was valid - if (_minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_minimizer->fitter()->Result().CovMatrixStatus() == 3) { break; } if (sIdx >= m_hessestrategy.Length() - 1) { @@ -1260,19 +1260,19 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, delete _nll; throw std::runtime_error("Keyboard interrupt while hesse calculating"); } - if ((_status != 0 || _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() != 3) && status == 0 && + if ((_status != 0 || _minimizer->fitter()->Result().CovMatrixStatus() != 3) && status == 0 && printLevel >= -1) { printCerr(TString::Format("Warning: %s hesse status is %d, covQual=%d", fitName.Data(), _status, - _minimizer->fitter()->GetMinimizer()->CovMatrixStatus()) + _minimizer->fitter()->Result().CovMatrixStatus()) .Data()); } - if (_status == 0 && _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_status == 0 && _minimizer->fitter()->Result().CovMatrixStatus() == 3) { // covariance is valid! break; } else if (_status == 0) { // set the statusHistory to the cov status, since that's more informative - statusHistory.back().second = _minimizer->fitter()->GetMinimizer()->CovMatrixStatus(); + statusHistory.back().second = _minimizer->fitter()->Result().CovMatrixStatus(); } if (sIdx >= m_hessestrategy.Length() - 1) {