44 :
RooAbsReal{
"RooEvaluatorWrapper",
"RooEvaluatorWrapper"},
48 _paramSet(
"paramSet",
"Set of parameters",
this),
54 setData(*
data,
false);
56 _paramSet.add(_evaluator->getParameters());
58 _paramSet.remove(*_paramSet.find(
item.first->GetName()));
64 _evaluator{
other._evaluator},
65 _topNode(
"topNode",
this,
other._topNode),
67 _paramSet(
"paramSet",
"Set of parameters",
this),
68 _rangeName{
other._rangeName},
70 _takeGlobalObservablesFromData{
other._takeGlobalObservablesFromData},
73 _paramSet.add(
other._paramSet);
76RooEvaluatorWrapper::~RooEvaluatorWrapper() =
default;
81 outputSet.add(_evaluator->getParameters());
83 outputSet.remove(*observables,
false,
true);
87 if (_data->getGlobalObservables() && _data->getGlobalObservables()->find(
item.first->GetName())) {
99 if (_takeGlobalObservablesFromData && _data->getGlobalObservables()) {
100 outputSet.replace(*_data->getGlobalObservables());
133 void gradient(
double *out)
const
136 std::fill(out, out + _params.size(), 0.0);
139 void hessian(
double *out)
const
142 std::fill(out, out + _params.size() * _params.size(), 0.0);
151 std::vector<std::string>
const &collectedFunctions() {
return _collectedFunctions; }
156 return _func(
_varBuffer.data(), _observables.data(), _xlArr.data());
167 using Func =
double (*)(
double *,
double const *,
double const *);
168 using Grad = void (*)(
double *,
double const *,
double const *,
double *);
169 using Hessian = void (*)(
double *,
double const *,
double const *,
double *);
172 std::string _funcName;
179 std::vector<double> _observables;
180 std::unordered_map<RooFit::Detail::DataKey, std::size_t>
_obsInfos;
181 std::vector<double> _xlArr;
182 std::vector<std::string> _collectedFunctions;
187void replaceAll(std::string &str,
const std::string &from,
const std::string &to)
193 str.replace(
start_pos, from.length(), to);
203 std::unordered_set<RooFit::Detail::DataKey> dependsOnData;
205 dependsOnData.insert(arg);
209 if (arg->getAttribute(
"__obs__")) {
210 dependsOnData.insert(arg);
213 if (
server->isValueServer(*arg)) {
214 if (dependsOnData.find(
server) != dependsOnData.end() && !arg->isReducerNode()) {
215 dependsOnData.insert(arg);
222 return dependsOnData;
244 std::unordered_set<RooFit::Detail::DataKey> dependsOnData;
255 ctx.
addResult(param,
"params[" + std::to_string(idx) +
"]");
266 auto print = [](std::string
const &
msg) {
oocoutI(
nullptr, Fitting) <<
msg << std::endl; };
271 gInterpreter->Declare(
"#include <RooFit/CodegenImpl.h>\n");
275 std::string
debugFileName =
"_codegen_" + _funcName +
".cxx";
276 errorMsg <<
"Function " << _funcName <<
" could not be compiled. See above for details. Full code dumped to file "
284 throw std::runtime_error(
errorMsg.str().c_str());
287 _func =
reinterpret_cast<Func
>(
gInterpreter->ProcessLine((_funcName +
";").c_str()));
289 _xlArr = ctx.
xlArr();
301 _observables.clear();
303 std::size_t
total = 0;
304 _observables.reserve(2 *
spans.size());
308 _observables.push_back(
total + 2 *
spans.size());
309 _observables.push_back(
item.second.size());
315 std::size_t
n =
item.second.size();
316 _observables.reserve(_observables.size() +
n);
317 for (std::size_t i = 0; i <
n; ++i) {
318 _observables.push_back(
item.second[i]);
324void RooFuncWrapper::createGradient()
327 std::string
gradName = _funcName +
"_grad_0";
331 gInterpreter->Declare(
"#include <Math/CladDerivator.h>\n");
337 " clad::gradient(" << _funcName <<
", \"params\");\n"
341 auto print = [](std::string
const &
msg) {
oocoutI(
nullptr, Fitting) <<
msg << std::endl; };
350 errorMsg <<
"Function could not be differentiated. See above for details.";
352 throw std::runtime_error(
errorMsg.str().c_str());
358 std::stringstream
ss;
360 ss <<
"static_cast<void (*)(double *, double const *, double const *, double *)>(" <<
gradName <<
");";
366 errorMsg <<
"Function could not be differentiated since ROOT was built without Clad support.";
368 throw std::runtime_error(
errorMsg.str().c_str());
372void RooFuncWrapper::createHessian()
375 std::string
hessianName = _funcName +
"_hessian_0";
376 std::string
requestName = _funcName +
"_hessian_req";
379 gInterpreter->Declare(
"#include <Math/CladDerivator.h>\n");
384 _params.size() == 1 ?
"\"params[0]\"" : (
"\"params[0:" + std::to_string(_params.size() - 1) +
"]\"");
387 " clad::hessian(" << _funcName <<
", " <<
paramsStr <<
");\n"
391 auto print = [](std::string
const &
msg) {
oocoutI(
nullptr, Fitting) <<
msg << std::endl; };
400 errorMsg <<
"Function could not be differentiated. See above for details.";
402 throw std::runtime_error(
errorMsg.str().c_str());
408 std::stringstream
ss;
410 ss <<
"static_cast<void (*)(double *, double const *, double const *, double *)>(" <<
hessianName <<
");";
416 errorMsg <<
"Function could not be differentiated since ROOT was built without Clad support.";
418 throw std::runtime_error(
errorMsg.str().c_str());
422void RooFuncWrapper::updateGradientVarBuffer()
const
425 return obj->isCategory() ? static_cast<RooAbsCategory *>(obj)->getCurrentIndex()
426 : static_cast<RooAbsReal *>(obj)->getVal();
431void RooFuncWrapper::writeDebugMacro(std::string
const &
filename)
const
437 for (std::string
const &
name : _collectedFunctions) {
442 std::unique_ptr<TInterpreterValue>
v =
gInterpreter->MakeInterpreterValue();
444 std::string s =
v->ToString();
445 for (
int i = 0; i < 2; ++i) {
446 s = s.erase(0, s.find(
"\n") + 1);
453 _params.size() == 1 ?
"\"params[0]\"" : (
"\"params[0:" + std::to_string(_params.size() - 1) +
"]\"");
455 outFile << R
"(//auto-generated test macro
456#include <RooFit/Detail/MathFuncs.h>
457#include <Math/CladDerivator.h>
464void gradient_request() {
466 << _funcName << R"(, "params");
469 << _funcName << ", " <<
paramsStr << R
"();
478 std::stringstream
decl;
479 decl <<
"std::vector<double> " <<
name <<
" = {";
480 for (std::size_t i = 0; i <
vec.size(); ++i) {
484 if (i <
vec.size() - 1)
497 outFile <<
"// clang-format off\n" << std::endl;
504 outFile <<
"// clang-format on\n" << std::endl;
507// To run as a ROOT macro
511 const std::size_t n = parametersVec.size();
513 std::vector<double> gradientVec(n);
515 auto func = [&](std::span<double> params) {
517 << _funcName << R"((params.data(), observablesVec.data(), auxConstantsVec.data());
519 auto grad = [&](std::span<double> params, std::span<double> out) {
521 << _funcName << R"(_grad_0(parametersVec.data(), observablesVec.data(), auxConstantsVec.data(),
525 grad(parametersVec, gradientVec);
527 auto numDiff = [&](int i) {
528 const double eps = 1e-6;
529 std::vector<double> p{parametersVec};
530 p[i] = parametersVec[i] - eps;
531 double funcValDown = func(p);
532 p[i] = parametersVec[i] + eps;
533 double funcValUp = func(p);
534 return (funcValUp - funcValDown) / (2 * eps);
537 for (std::size_t i = 0; i < parametersVec.size(); ++i) {
538 std::cout << i << ":" << std::endl;
539 std::cout << " numr : " << numDiff(i) << std::endl;
540 std::cout << " clad : " << gradientVec[i] << std::endl;
546 auto hess = [&](std::span<double> params, std::span<double> out) {
548 << _funcName << R"(_hessian_0(params.data(), observablesVec.data(), auxConstantsVec.data(), out.data());
551 std::vector<double> hessianVec(n * n);
552 hess(parametersVec, hessianVec);
554 // ---------- Numerical Hessian ----------
555 // Uses central differences:
556 // diag: (f(x+ei)-2f(x)+f(x-ei))/eps^2
557 // offdiag: (f(++ ) - f(+-) - f(-+) + f(--)) / (4 eps^2)
558 auto numHess = [&](std::size_t i, std::size_t j) {
559 const double eps = 1e-5; // often needs to be a bit larger than grad eps
560 std::vector<double> p(parametersVec.begin(), parametersVec.end());
563 const double f0 = func(p);
565 p[i] = parametersVec[i] + eps;
566 const double fUp = func(p);
568 p[i] = parametersVec[i] - eps;
569 const double fDown = func(p);
571 return (fUp - 2.0 * f0 + fDown) / (eps * eps);
573 // f(x_i + eps, x_j + eps)
574 p[i] = parametersVec[i] + eps;
575 p[j] = parametersVec[j] + eps;
576 const double fPP = func(p);
578 // f(x_i + eps, x_j - eps)
579 p[i] = parametersVec[i] + eps;
580 p[j] = parametersVec[j] - eps;
581 const double fPM = func(p);
583 // f(x_i - eps, x_j + eps)
584 p[i] = parametersVec[i] - eps;
585 p[j] = parametersVec[j] + eps;
586 const double fMP = func(p);
588 // f(x_i - eps, x_j - eps)
589 p[i] = parametersVec[i] - eps;
590 p[j] = parametersVec[j] - eps;
591 const double fMM = func(p);
593 return (fPP - fPM - fMP + fMM) / (4.0 * eps * eps);
597 // Compute full numerical Hessian
598 std::vector<double> numHessianVec(n * n);
599 for (std::size_t i = 0; i < n; ++i) {
600 for (std::size_t j = 0; j < n; ++j) {
601 numHessianVec[i + n * j] = numHess(i, j); // keep same layout as your print
605 // ---------- Compare & print ----------
606 std::cout << "Hessian comparison (clad vs numeric vs diff):\n\n";
608 for (std::size_t i = 0; i < n; ++i) {
609 for (std::size_t j = 0; j < n; ++j) {
610 const std::size_t idx = i + n * j; // same indexing you used
611 const double cladH = hessianVec[idx];
612 const double numH = numHessianVec[idx];
613 const double diff = cladH - numH;
615 std::cout << "[" << i << "," << j << "] "
616 << "clad=" << cladH << " num=" << numH << " diff=" << diff << "\n";
620 std::cout << "\nRaw Clad Hessian matrix:\n";
621 for (std::size_t i = 0; i < n; ++i) {
622 for (std::size_t j = 0; j < n; ++j) {
623 std::cout << hessianVec[i + n * j] << " ";
628 std::cout << "\nRaw Numerical Hessian matrix:\n";
629 for (std::size_t i = 0; i < n; ++i) {
630 for (std::size_t j = 0; j < n; ++j) {
631 std::cout << numHessianVec[i + n * j] << " ";
640double RooEvaluatorWrapper::evaluate()
const
649 :
RooFit::EvalContext::OffsetMode::WithOffset);
651 return _evaluator->run()[0];
660 constexpr auto errMsg =
"Error in RooAbsReal::setData(): only resetting with same-structured data is supported.";
666 std::stack<std::vector<double>>{}.swap(_vectorBuffers);
667 const bool isChi2 = _topNode->getAttribute(
"Chi2EvaluationActive");
671 _takeGlobalObservablesFromData, _vectorBuffers);
674 throw std::runtime_error(
errMsg);
677 const char *
name =
item.first->GetName();
678 _evaluator->setInput(
name,
item.second,
false);
679 if (_paramSet.find(
name)) {
681 throw std::runtime_error(
errMsg);
690void RooEvaluatorWrapper::createFuncWrapper()
694 this->getParameters(_data ? _data->get() : nullptr,
paramSet,
false);
696 const bool isChi2 = _topNode->getAttribute(
"Chi2EvaluationActive");
702void RooEvaluatorWrapper::generateGradient()
710void RooEvaluatorWrapper::generateHessian()
718void RooEvaluatorWrapper::setUseGeneratedFunctionCode(
bool flag)
725void RooEvaluatorWrapper::gradient(
double *out)
const
730void RooEvaluatorWrapper::hessian(
double *out)
const
735bool RooEvaluatorWrapper::hasGradient()
const
740bool RooEvaluatorWrapper::hasHessian()
const
745void RooEvaluatorWrapper::writeDebugMacro(std::string
const &
filename)
const
753 return _evaluator->setOperModes(
opMode);
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
static unsigned int total
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
const_iterator begin() const
const_iterator end() const
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Abstract base class for binned and unbinned datasets.
Abstract interface for all probability density functions.
Abstract base class for objects that represent a real value and implements functionality common to al...
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
A class to maintain the context for squashing of RooFit models into code.
void addResult(RooAbsArg const *key, std::string const &value)
A function to save an expression that includes/depends on the result of the input node.
void addVecObs(const char *key, int idx)
Since the squashed code represents all observables as a single flattened array, it is important to ke...
std::string buildFunction(RooAbsArg const &arg, std::unordered_set< RooFit::Detail::DataKey > const &dependsOnData={})
Assemble and return the final code with the return expression and global statements.
std::vector< std::string > const & collectedFunctions()
std::string const & collectedCode()
std::vector< double > const & xlArr()
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
void replaceAll(std::string &inOut, std::string_view what, std::string_view with)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
void getSortedComputationGraph(RooAbsArg const &func, RooArgSet &out)
void evaluate(typename Architecture_t::Tensor_t &A, EActivationFunction f)
Apply the given activation function to each value in the given tensor A.