From 40f65c4b2c8a34e87ef1b2859d7e13e945dc7ad5 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 27 Aug 2026 15:28:24 +0000 Subject: [PATCH 1/5] [RF] Add JIT-free AST evaluator for RooFormula MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add a small, self-contained expression compiler and evaluator to RooFit (RooFormulaParser + RooExprEvaluator) that handles the realistic subset of RooFormula expressions without any use of the interpreter/JIT. This is the main building block for fixing cling out-of-memory problems with many expression-based NormFactors (see #21052; not fully closed yet, since the codegen path still JIT-compiles on demand). The contract is strict: either the compiled program evaluates bitwise identically to what TFormula/cling computes for the same processed formula string, or the parser refuses and RooFormula silently falls back to the TFormula backend, so every expression that worked before keeps working, with identical values and identical error messages for invalid formulas. In particular: - Operator precedence matches C++; `^`/`**` reproduce TFormula's HandleExponentiation rewrite (right-associative exponentiation binding tighter than `*` and unary minus, `expr^2` -> TMath::Sq). - cling's int-vs-double expression typing is tracked: truncating integer division (`1/2`, `(x>0)/2`) and min/max with mixed int/double arguments (invalid in cling) fall back. - Function calls evaluate exactly what the JIT'd code called for each spelling: bare/std:: spellings call libm/std, TMath:: spellings call the TMath functions (which are not always identical, e.g. TMath::Erf, TMath::ATan2, TMath::Min). - Comparisons yield exactly 0.0/1.0; `&&`/`||` do not short-circuit and `?:` evaluates both branches (keeping the scalar semantics consistent with a future vectorized path). Identical formula strings share one immutable instruction vector through a mutex-protected process-wide registry (mirroring what TFormula's gClingFunctions cache does for JIT'd code); evaluation is const, lock-free and touches no globals, so it is safe under RooFit's concurrent evaluation. The environment variable ROOFIT_FORMULA_BACKEND overrides the backend choice: `tformula` never uses the new evaluator, `ast` fails loudly instead of falling back (for testing), unset means AST with silent fallback. The first fallback in a process is reported once at INFO level, so it is not completely invisible; the individual fallbacks stay on the InputArguments debug stream. Since the codegen path calls the cling-JIT-compiled TFormula function by name (getUniqueFuncName()), RooFormula::getTFormula() now lazily creates a TFormula on demand when the formula is evaluated JIT-free, keeping codegen and AD working exactly as before. This stopgap goes away when codegen emits C++ from the parsed expression directly. On the real-world formula corpus (all RooFormulaVar/RooGenericPdf strings in the RooFit tests and tutorials, RooFit-generated formulas, and the RooFit-realistic subset of the TFormula parsing tests), 152 of 156 entries take the JIT-free path (97%) and evaluate bitwise identically to TFormula on random inputs; the remainder is the ROOT::Math::*_pdf family, which falls back. 🤖 Done with the help of AI --- roofit/roofitcore/CMakeLists.txt | 6 +- roofit/roofitcore/inc/RooFormulaVar.h | 2 + roofit/roofitcore/inc/RooGenericPdf.h | 2 + roofit/roofitcore/src/RooExprEvaluator.cxx | 345 +++++++ roofit/roofitcore/src/RooExprEvaluator.h | 132 +++ roofit/roofitcore/src/RooFormulaParser.cxx | 808 +++++++++++++++++ roofit/roofitcore/src/RooFormulaParser.h | 44 + roofit/roofitcore/src/RooFormulaUtils.cxx | 116 +++ roofit/roofitcore/src/RooFormulaUtils.h | 8 + roofit/roofitcore/src/RooFormulaVar.cxx | 23 +- roofit/roofitcore/src/RooGenericPdf.cxx | 15 +- roofit/roofitcore/test/CMakeLists.txt | 1 + .../test/testRooFormulaEvaluator.cxx | 848 ++++++++++++++++++ 13 files changed, 2347 insertions(+), 3 deletions(-) create mode 100644 roofit/roofitcore/src/RooExprEvaluator.cxx create mode 100644 roofit/roofitcore/src/RooExprEvaluator.h create mode 100644 roofit/roofitcore/src/RooFormulaParser.cxx create mode 100644 roofit/roofitcore/src/RooFormulaParser.h create mode 100644 roofit/roofitcore/test/testRooFormulaEvaluator.cxx diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 6e7429abe6164..a0a6e556b3dac 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -343,6 +343,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/RooErrorVar.cxx src/RooEvaluatorWrapper.cxx src/RooExpensiveObjectCache.cxx + src/RooExprEvaluator.cxx src/RooExtendPdf.cxx src/RooExtendedBinding.cxx src/RooExtendedTerm.cxx @@ -356,6 +357,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/RooFitLegacy/RooCatTypeLegacy.cxx src/RooFitResult.cxx src/RooFoamGenerator.cxx + src/RooFormulaParser.cxx src/RooFormulaUtils.cxx src/RooFormulaVar.cxx src/RooFracRemainder.cxx @@ -540,8 +542,10 @@ if(NOT CMAKE_VERSION VERSION_LESS "3.23.0") inc/LinkDef.h inc/RooFit/UniqueId.h # Was not being included in ROOT_STL_PACKAGE call src/RooConvIntegrandBinding.h - src/RooFormulaUtils.h + src/RooExprEvaluator.h src/RooFormulaEvaluator.h + src/RooFormulaParser.h + src/RooFormulaUtils.h src/RooTFormulaEvaluator.h src/RooFitLegacy/RooAbsCategoryLegacyIterator.h src/RooMCIntegrator.h diff --git a/roofit/roofitcore/inc/RooFormulaVar.h b/roofit/roofitcore/inc/RooFormulaVar.h index 047c6455bca04..b7cf0232090a3 100644 --- a/roofit/roofitcore/inc/RooFormulaVar.h +++ b/roofit/roofitcore/inc/RooFormulaVar.h @@ -29,6 +29,7 @@ class RooArgSet ; class RooFormulaEvaluator; class RooAbsRealLValue; +class TFormula; class RooFormulaVar : public RooAbsReal { public: @@ -96,6 +97,7 @@ class RooFormulaVar : public RooAbsReal { RooListProxy _actualVars; ///< Actual parameters used by formula engine mutable std::unique_ptr _evaluator; /// _tFormulaForCodegen; /// _evaluator; /// _tFormulaForCodegen; ///> _binnings; ///< User-defined binnings, keyed by the observable's index diff --git a/roofit/roofitcore/src/RooExprEvaluator.cxx b/roofit/roofitcore/src/RooExprEvaluator.cxx new file mode 100644 index 0000000000000..51695d0bc5554 --- /dev/null +++ b/roofit/roofitcore/src/RooExprEvaluator.cxx @@ -0,0 +1,345 @@ +/// \cond ROOFIT_INTERNAL + +/* + * Project: RooFit + * + * Copyright (c) 2026, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +#include "RooExprEvaluator.h" + +#include "TMath.h" + +#include +#include +#include +#include + +namespace RooFormulaFunctions { + +namespace { + +// Entry construction helpers (C++17 has no designated initializers). +Entry F0(const char *name, double (*fn)()) +{ + Entry e; + e.name = name; + e.arity = 0; + e.fn0 = fn; + return e; +} + +Entry F1(const char *name, double (*fn)(double), TypeRule rule = TypeRule::Double) +{ + Entry e; + e.name = name; + e.arity = 1; + e.rule = rule; + e.fn1 = fn; + return e; +} + +Entry F2(const char *name, double (*fn)(double, double), TypeRule rule = TypeRule::Double) +{ + Entry e; + e.name = name; + e.arity = 2; + e.rule = rule; + e.fn2 = fn; + return e; +} + +Entry F3(const char *name, double (*fn)(double, double, double)) +{ + Entry e; + e.name = name; + e.arity = 3; + e.fn3 = fn; + return e; +} + +Entry F4(const char *name, double (*fn)(double, double, double, double)) +{ + Entry e; + e.name = name; + e.arity = 4; + e.fn4 = fn; + return e; +} + +// The allow-list. For every accepted spelling, the function pointer calls +// exactly what the cling-JIT-compiled TFormula code called for that spelling: +// bare and std:: spellings resolve to the libm/std functions, TMath:: +// spellings call the TMath functions (whose implementations are not always +// identical to libm, e.g. TMath::Erf goes through ROOT::Math::erf and +// TMath::ATan2 special-cases x == 0). This makes the two evaluation backends +// agree bitwise. +// +// Bare spellings work in TFormula because cling resolves them against the +// global namespace (libm) and `using namespace std`; the only genuine +// TFormula shortcuts are `sign` -> TMath::Sign and `sq` -> TMath::Sq +// (see TFormula::FillDefaults). All spellings below were cross-checked +// against what TFormula accepts today; do not add a spelling TFormula +// would reject (e.g. there is no std::sq or std::sign). +std::vector makeTable() +{ + auto castInt = +[](double x) { return static_cast(static_cast(x)); }; + auto square = +[](double x) { return x * x; }; + // std::min/std::max compare exactly like this; the asymmetric NaN behavior + // (min(NaN, 1) = NaN but min(1, NaN) = 1) must be reproduced, so no fmin/fmax. + auto stdMin = +[](double a, double b) { return b < a ? b : a; }; + auto stdMax = +[](double a, double b) { return a < b ? b : a; }; + // TMath::Min/Max compare differently from std::min/max (visible with NaNs). + auto tmathMin = +[](double a, double b) { return TMath::Min(a, b); }; + auto tmathMax = +[](double a, double b) { return TMath::Max(a, b); }; + // Both the `sign` shortcut and TMath::Sign resolve to TMath::Sign, which is + // std::copysign for double arguments. (For an integer first argument cling + // picks the TMath::Sign template, whose result is numerically identical.) + auto sign = +[](double a, double b) { return TMath::Sign(a, b); }; + auto signBit = +[](double x) { return std::signbit(x) ? 1.0 : 0.0; }; + + return { + // clang-format off + // one-argument functions, libm/std spellings + F1("sqrt", +[](double x) { return std::sqrt(x); }), + F1("std::sqrt", +[](double x) { return std::sqrt(x); }), + F1("exp", +[](double x) { return std::exp(x); }), + F1("std::exp", +[](double x) { return std::exp(x); }), + F1("log", +[](double x) { return std::log(x); }), + F1("std::log", +[](double x) { return std::log(x); }), + F1("log10", +[](double x) { return std::log10(x); }), + F1("std::log10", +[](double x) { return std::log10(x); }), + F1("sin", +[](double x) { return std::sin(x); }), + F1("std::sin", +[](double x) { return std::sin(x); }), + F1("cos", +[](double x) { return std::cos(x); }), + F1("std::cos", +[](double x) { return std::cos(x); }), + F1("tan", +[](double x) { return std::tan(x); }), + F1("std::tan", +[](double x) { return std::tan(x); }), + F1("asin", +[](double x) { return std::asin(x); }), + F1("std::asin", +[](double x) { return std::asin(x); }), + F1("acos", +[](double x) { return std::acos(x); }), + F1("std::acos", +[](double x) { return std::acos(x); }), + F1("atan", +[](double x) { return std::atan(x); }), + F1("std::atan", +[](double x) { return std::atan(x); }), + F1("sinh", +[](double x) { return std::sinh(x); }), + F1("std::sinh", +[](double x) { return std::sinh(x); }), + F1("cosh", +[](double x) { return std::cosh(x); }), + F1("std::cosh", +[](double x) { return std::cosh(x); }), + F1("tanh", +[](double x) { return std::tanh(x); }), + F1("std::tanh", +[](double x) { return std::tanh(x); }), + F1("asinh", +[](double x) { return std::asinh(x); }), + F1("std::asinh", +[](double x) { return std::asinh(x); }), + F1("acosh", +[](double x) { return std::acosh(x); }), + F1("std::acosh", +[](double x) { return std::acosh(x); }), + F1("atanh", +[](double x) { return std::atanh(x); }), + F1("std::atanh", +[](double x) { return std::atanh(x); }), + F1("floor", +[](double x) { return std::floor(x); }), + F1("std::floor", +[](double x) { return std::floor(x); }), + F1("ceil", +[](double x) { return std::ceil(x); }), + F1("std::ceil", +[](double x) { return std::ceil(x); }), + F1("erf", +[](double x) { return std::erf(x); }), + F1("std::erf", +[](double x) { return std::erf(x); }), + F1("erfc", +[](double x) { return std::erfc(x); }), + F1("std::erfc", +[](double x) { return std::erfc(x); }), + F1("tgamma", +[](double x) { return std::tgamma(x); }), + F1("std::tgamma", +[](double x) { return std::tgamma(x); }), + F1("lgamma", +[](double x) { return std::lgamma(x); }), + F1("std::lgamma", +[](double x) { return std::lgamma(x); }), + // abs is integer-preserving in C++ (::abs(int), std::abs(int)); fabs is not + F1("abs", +[](double x) { return std::fabs(x); }, TypeRule::SameAsFirstArg), + F1("std::abs", +[](double x) { return std::fabs(x); }, TypeRule::SameAsFirstArg), + F1("fabs", +[](double x) { return std::fabs(x); }), + F1("std::fabs", +[](double x) { return std::fabs(x); }), + F1("int", castInt, TypeRule::Int), // C++ functional cast: truncation towards zero + F1("sq", square), // TFormula shortcut for TMath::Sq(Double_t) + // one-argument functions, TMath spellings + F1("TMath::Sqrt", +[](double x) { return TMath::Sqrt(x); }), + F1("TMath::Exp", +[](double x) { return TMath::Exp(x); }), + F1("TMath::Log", +[](double x) { return TMath::Log(x); }), + F1("TMath::Log10", +[](double x) { return TMath::Log10(x); }), + F1("TMath::Sin", +[](double x) { return TMath::Sin(x); }), + F1("TMath::Cos", +[](double x) { return TMath::Cos(x); }), + F1("TMath::Tan", +[](double x) { return TMath::Tan(x); }), + F1("TMath::ASin", +[](double x) { return TMath::ASin(x); }), + F1("TMath::ACos", +[](double x) { return TMath::ACos(x); }), + F1("TMath::ATan", +[](double x) { return TMath::ATan(x); }), + F1("TMath::SinH", +[](double x) { return TMath::SinH(x); }), + F1("TMath::CosH", +[](double x) { return TMath::CosH(x); }), + F1("TMath::TanH", +[](double x) { return TMath::TanH(x); }), + F1("TMath::ASinH", +[](double x) { return TMath::ASinH(x); }), + F1("TMath::ACosH", +[](double x) { return TMath::ACosH(x); }), + F1("TMath::ATanH", +[](double x) { return TMath::ATanH(x); }), + F1("TMath::Floor", +[](double x) { return TMath::Floor(x); }), + F1("TMath::Ceil", +[](double x) { return TMath::Ceil(x); }), + F1("TMath::Erf", +[](double x) { return TMath::Erf(x); }), + F1("TMath::Erfc", +[](double x) { return TMath::Erfc(x); }), + F1("TMath::Abs", +[](double x) { return TMath::Abs(x); }, TypeRule::SameAsFirstArg), + F1("TMath::Sq", square), + F1("TMath::SignBit", signBit, TypeRule::Int), + // two-argument functions + F2("pow", +[](double a, double b) { return std::pow(a, b); }), + F2("std::pow", +[](double a, double b) { return std::pow(a, b); }), + F2("TMath::Power", +[](double a, double b) { return TMath::Power(a, b); }), + F2("atan2", +[](double a, double b) { return std::atan2(a, b); }), + F2("std::atan2", +[](double a, double b) { return std::atan2(a, b); }), + F2("TMath::ATan2", +[](double a, double b) { return TMath::ATan2(a, b); }), + F2("fmod", +[](double a, double b) { return std::fmod(a, b); }), + F2("std::fmod", +[](double a, double b) { return std::fmod(a, b); }), + F2("min", stdMin, TypeRule::MinMax), + F2("std::min", stdMin, TypeRule::MinMax), + F2("TMath::Min", tmathMin, TypeRule::MinMax), + F2("max", stdMax, TypeRule::MinMax), + F2("std::max", stdMax, TypeRule::MinMax), + F2("TMath::Max", tmathMax, TypeRule::MinMax), + F2("sign", sign, TypeRule::SameAsFirstArg), // TFormula shortcut for TMath::Sign + F2("TMath::Sign", sign, TypeRule::SameAsFirstArg), + // zero-argument constants (folded to Op::Const at parse time) + F0("TMath::Pi", +[]() { return TMath::Pi(); }), + F0("TMath::TwoPi", +[]() { return TMath::TwoPi(); }), + F0("TMath::PiOver2",+[]() { return TMath::PiOver2(); }), + F0("TMath::E", +[]() { return TMath::E(); }), + // TMath::Gaus with its default arguments mean=0, sigma=1, norm=false + F1("TMath::Gaus", +[](double x) { return TMath::Gaus(x); }), + F2("TMath::Gaus", +[](double x, double m) { return TMath::Gaus(x, m); }), + F3("TMath::Gaus", +[](double x, double m, double s) { return TMath::Gaus(x, m, s); }), + F4("TMath::Gaus", +[](double x, double m, double s, double n) { return TMath::Gaus(x, m, s, n != 0.0); }), + // clang-format on + }; +} + +std::vector const &theTable() +{ + static const std::vector t = makeTable(); + return t; +} + +} // namespace + +Entry const *table() +{ + return theTable().data(); +} + +std::size_t tableSize() +{ + return theTable().size(); +} + +Entry const *find(std::string const &name, unsigned int nArgs, std::uint32_t &index) +{ + auto const &tab = theTable(); + for (std::size_t i = 0; i < tab.size(); ++i) { + if (tab[i].arity == nArgs && name == tab[i].name) { + index = static_cast(i); + return &tab[i]; + } + } + return nullptr; +} + +} // namespace RooFormulaFunctions + +//////////////////////////////////////////////////////////////////////////////// +/// Interpret the instruction sequence with a fixed-size value stack. Stack +/// depth was bounded at parse time and is re-checked once up front, so no +/// per-instruction bounds checks are needed. +double RooExprEvaluator::eval(const double *vars) const +{ + // The depth bound is established at parse time; refuse to evaluate rather + // than trust it, so a future accounting bug cannot overflow the stack. + if (_program->stackDepth > kMaxStackDepth) { + throw std::runtime_error("RooExprEvaluator: expression program exceeds the maximum stack depth"); + } + + double stack[kMaxStackDepth]; + std::size_t sp = 0; + auto const *funcs = RooFormulaFunctions::table(); + + for (Instr const &ins : _program->code) { + switch (ins.op) { + case Op::Const: stack[sp++] = ins.konst; break; + case Op::Var: stack[sp++] = vars[ins.arg]; break; + case Op::Add: + --sp; + stack[sp - 1] = stack[sp - 1] + stack[sp]; + break; + case Op::Sub: + --sp; + stack[sp - 1] = stack[sp - 1] - stack[sp]; + break; + case Op::Mul: + --sp; + stack[sp - 1] = stack[sp - 1] * stack[sp]; + break; + case Op::Div: + --sp; + stack[sp - 1] = stack[sp - 1] / stack[sp]; + break; + case Op::Neg: stack[sp - 1] = -stack[sp - 1]; break; + case Op::Not: stack[sp - 1] = stack[sp - 1] == 0.0 ? 1.0 : 0.0; break; + case Op::LT: + --sp; + stack[sp - 1] = stack[sp - 1] < stack[sp] ? 1.0 : 0.0; + break; + case Op::LE: + --sp; + stack[sp - 1] = stack[sp - 1] <= stack[sp] ? 1.0 : 0.0; + break; + case Op::GT: + --sp; + stack[sp - 1] = stack[sp - 1] > stack[sp] ? 1.0 : 0.0; + break; + case Op::GE: + --sp; + stack[sp - 1] = stack[sp - 1] >= stack[sp] ? 1.0 : 0.0; + break; + case Op::EQ: + --sp; + stack[sp - 1] = stack[sp - 1] == stack[sp] ? 1.0 : 0.0; + break; + case Op::NE: + --sp; + stack[sp - 1] = stack[sp - 1] != stack[sp] ? 1.0 : 0.0; + break; + case Op::And: + --sp; + stack[sp - 1] = (stack[sp - 1] != 0.0 && stack[sp] != 0.0) ? 1.0 : 0.0; + break; + case Op::Or: + --sp; + stack[sp - 1] = (stack[sp - 1] != 0.0 || stack[sp] != 0.0) ? 1.0 : 0.0; + break; + case Op::Select: + sp -= 2; + stack[sp - 1] = stack[sp - 1] != 0.0 ? stack[sp] : stack[sp + 1]; + break; + case Op::Pow: + --sp; + stack[sp - 1] = std::pow(stack[sp - 1], stack[sp]); + break; + case Op::Sq: stack[sp - 1] *= stack[sp - 1]; break; + case Op::IntNorm: stack[sp - 1] += 0.0; break; + case Op::Call1: stack[sp - 1] = funcs[ins.arg].fn1(stack[sp - 1]); break; + case Op::Call2: + --sp; + stack[sp - 1] = funcs[ins.arg].fn2(stack[sp - 1], stack[sp]); + break; + case Op::Call3: + sp -= 2; + stack[sp - 1] = funcs[ins.arg].fn3(stack[sp - 1], stack[sp], stack[sp + 1]); + break; + case Op::Call4: + sp -= 3; + stack[sp - 1] = funcs[ins.arg].fn4(stack[sp - 1], stack[sp], stack[sp + 1], stack[sp + 2]); + break; + } + } + + return stack[0]; +} + +/// \endcond diff --git a/roofit/roofitcore/src/RooExprEvaluator.h b/roofit/roofitcore/src/RooExprEvaluator.h new file mode 100644 index 0000000000000..45d743199940b --- /dev/null +++ b/roofit/roofitcore/src/RooExprEvaluator.h @@ -0,0 +1,132 @@ +/// \cond ROOFIT_INTERNAL + +/* + * Project: RooFit + * + * Copyright (c) 2026, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +#ifndef ROO_EXPR_EVALUATOR +#define ROO_EXPR_EVALUATOR + +#include "RooFormulaEvaluator.h" + +#include +#include +#include +#include + +/// Functions that can be called from a formula on the JIT-free evaluation +/// path. Each table entry corresponds to one accepted spelling (e.g. `sin`, +/// `std::sin`, and `TMath::Sin` are three entries), and its function pointer +/// reproduces exactly what the cling-JIT-compiled code called for that +/// spelling, so that both backends agree bitwise. +namespace RooFormulaFunctions { + +/// How the C++ result type of a call depends on the argument types. The parser +/// tracks int-ness of subexpressions to reproduce cling's expression typing +/// (in particular to detect integer division, which is not supported). +enum class TypeRule : std::uint8_t { + Double, ///< result is always double + SameAsFirstArg, ///< result type equals the type of the first argument (abs, sign) + Int, ///< result is an integer type (`int(x)` cast, TMath::SignBit) + MinMax ///< int if both args are int; mixed int/double does not compile in cling +}; + +struct Entry { + const char *name = nullptr; ///< accepted spelling in the formula + std::uint8_t arity = 0; + TypeRule rule = TypeRule::Double; + double (*fn0)() = nullptr; + double (*fn1)(double) = nullptr; + double (*fn2)(double, double) = nullptr; + double (*fn3)(double, double, double) = nullptr; + double (*fn4)(double, double, double, double) = nullptr; +}; + +Entry const *table(); +std::size_t tableSize(); + +/// Find the entry for the given spelling and argument count, or return +/// nullptr if there is none. The index output is the position in table(). +Entry const *find(std::string const &name, unsigned int nArgs, std::uint32_t &index); + +} // namespace RooFormulaFunctions + +/// RooFormulaEvaluator implementation that interprets a compiled postfix +/// instruction sequence, without any use of the interpreter/JIT. Instances are +/// created via RooFormulaParser::compile(), which returns a shared immutable +/// Program: identical formula strings share one instruction vector. +/// +/// eval() is const, touches no globals and no mutable state, so it is safe to +/// call concurrently from multiple threads. +class RooExprEvaluator final : public RooFormulaEvaluator { +public: + enum class Op : std::uint8_t { + Const, ///< push konst + Var, ///< push vars[arg] + Add, ///< a + b + Sub, ///< a - b + Mul, ///< a * b + Div, ///< a / b + Neg, ///< -a + Not, ///< !a (exactly 0.0 or 1.0) + LT, ///< a < b (exactly 0.0 or 1.0, likewise below) + LE, ///< a <= b + GT, ///< a > b + GE, ///< a >= b + EQ, ///< a == b + NE, ///< a != b + And, ///< a && b (no short-circuit: both operands are always evaluated) + Or, ///< a || b (no short-circuit) + Select, ///< c ? a : b (both branches are always evaluated) + Pow, ///< std::pow(a, b), from the `^`/`**` operator or pow() + Sq, ///< a * a, from TFormula's `expr^2` -> TMath::Sq(expr) rewrite + IntNorm, ///< a + 0.0: maps -0.0 to +0.0 where cling would have used integer arithmetic + Call1, ///< RooFormulaFunctions::table()[arg].fn1 + Call2, ///< RooFormulaFunctions::table()[arg].fn2 + Call3, ///< RooFormulaFunctions::table()[arg].fn3 + Call4 ///< RooFormulaFunctions::table()[arg].fn4 + }; + + struct Instr { + Op op = Op::Const; + std::uint32_t arg = 0; + double konst = 0.0; + }; + + /// A compiled formula: immutable after construction and shared between all + /// RooFormula instances with the same processed formula string. + struct Program { + std::vector code; + std::vector usedVars; ///< usedVars[i] is true if `x[i]` appears in the formula + std::string formula; ///< the processed formula string this was compiled from + unsigned int stackDepth = 0; + }; + + /// Maximum evaluation stack depth (checked at compile time in the parser). + static constexpr unsigned int kMaxStackDepth = 256; + + RooExprEvaluator(std::shared_ptr program) : _program{std::move(program)} {} + + double eval(const double *vars) const override; + + std::unique_ptr clone() const override { return std::make_unique(_program); } + + /// Whether `x[i]` appears in the formula, as recorded while parsing. + bool usesVariable(unsigned int i) const { return i < _program->usedVars.size() && _program->usedVars[i]; } + + /// The processed formula string this program was compiled from. + std::string processedFormula() const { return _program->formula; } + +private: + std::shared_ptr _program; +}; + +#endif + +/// \endcond diff --git a/roofit/roofitcore/src/RooFormulaParser.cxx b/roofit/roofitcore/src/RooFormulaParser.cxx new file mode 100644 index 0000000000000..9635bab231320 --- /dev/null +++ b/roofit/roofitcore/src/RooFormulaParser.cxx @@ -0,0 +1,808 @@ +/// \cond ROOFIT_INTERNAL + +/* + * Project: RooFit + * + * Copyright (c) 2026, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +/** + * Recursive-descent parser that compiles a processed RooFormula expression + * string into the postfix instruction sequence of RooExprEvaluator. + * + * The contract is: either the compiled program evaluates *bitwise identically* + * to what TFormula/cling computes for the same string, or compile() fails and + * the caller falls back to the TFormula backend. Consequently: + * + * - The tokenizer owns all character-level concerns (numbers, multi-character + * operators, `x[i]` variables, `::`-qualified names); the parser only deals + * in tokens. + * - Operator precedence matches C++, which is what cling compiled. The one + * deliberate dialect difference is `^` (and `**`), which TFormula rewrites + * to pow()/TMath::Sq() *before* cling sees the string: it is + * right-associative exponentiation binding tighter than unary minus, whose + * right-hand side may carry one leading sign (see + * TFormula::HandleExponentiation). + * - cling's expression typing is tracked (integer literals and bool-valued + * operators are `int` in C++): integer division like `1/2` or `(x>0)/2` + * truncates in cling, so such expressions are not supported here and fall + * back. Similarly min/max with mixed int/double arguments does not compile + * in cling at all and is rejected. + * - `%` on doubles does not compile in cling, so TFormula formulas using it + * are invalid today; it is not part of this grammar either. + * - `&&` and `||` do not short-circuit and `?:` evaluates both branches, so + * that the scalar and a future vectorized path behave identically. The + * selected/combined values are unchanged; this is only observable if + * floating-point exceptions are trapped. + */ + +#include "RooFormulaParser.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using Op = RooExprEvaluator::Op; +using Instr = RooExprEvaluator::Instr; +using Program = RooExprEvaluator::Program; + +enum class Tok : std::uint8_t { + Number, + Var, + Ident, + LParen, + RParen, + Comma, + Plus, + Minus, + Star, + Slash, + Caret, + Lt, + Le, + Gt, + Ge, + Eq, + Ne, + AndAnd, + OrOr, + Not, + Question, + Colon, + End +}; + +struct Token { + Tok kind = Tok::End; + double value = 0.0; ///< Tok::Number + bool isInt = false; ///< Tok::Number: literal has integer type in C++ + std::uint32_t varIndex = 0; ///< Tok::Var + std::string_view text; ///< Tok::Ident / Tok::Number: raw spelling +}; + +bool isIdentStart(char c) +{ + return std::isalpha(static_cast(c)) || c == '_'; +} + +bool isIdentChar(char c) +{ + return std::isalnum(static_cast(c)) || c == '_'; +} + +class Tokenizer { +public: + Tokenizer(std::string const &s, std::string &error) : _s{s}, _error{error} {} + + bool run(std::vector &out) + { + const std::size_t n = _s.size(); + std::size_t i = 0; + while (i < n) { + const char c = _s[i]; + if (std::isspace(static_cast(c))) { + ++i; + continue; + } + if (std::isdigit(static_cast(c)) || + (c == '.' && i + 1 < n && std::isdigit(static_cast(_s[i + 1])))) { + if (!lexNumber(i, out)) + return false; + continue; + } + if (isIdentStart(c)) { + if (!lexIdentOrVar(i, out)) + return false; + continue; + } + if (!lexOperator(i, out)) + return false; + } + out.push_back(Token{}); + out.back().kind = Tok::End; + return true; + } + +private: + bool fail(std::string msg) + { + _error = std::move(msg); + return false; + } + + /// Lex one numeric literal, reproducing the C++ literal type: a literal is + /// `int` unless it has a decimal point or a (well-formed) exponent. Octal + /// (`010`) and hex (`0x64`) integer literals are handled by strtoll with + /// base 0, exactly as cling reads them. Literal values are parsed with + /// strtod/strtoll, never by hand, so they are correctly rounded like the + /// compiler's. + bool lexNumber(std::size_t &i, std::vector &out) + { + const std::size_t n = _s.size(); + const char *begin = _s.c_str() + i; + char *end = nullptr; + + bool isFloat = false; + const bool isHex = _s[i] == '0' && i + 1 < n && (_s[i + 1] == 'x' || _s[i + 1] == 'X'); + if (!isHex) { + std::size_t j = i; + while (j < n && std::isdigit(static_cast(_s[j]))) + ++j; + if (j < n && _s[j] == '.') { + isFloat = true; + } else if (j < n && (_s[j] == 'e' || _s[j] == 'E')) { + std::size_t k = j + 1; + if (k < n && (_s[k] == '+' || _s[k] == '-')) + ++k; + if (k < n && std::isdigit(static_cast(_s[k]))) + isFloat = true; + // else: something like "1e" -- fails the trailing-character check below + } + } + + Token tok; + tok.kind = Tok::Number; + errno = 0; + if (isFloat) { + tok.value = std::strtod(begin, &end); + tok.isInt = false; + } else { + const long long v = std::strtoll(begin, &end, 0); + if (errno == ERANGE) + return fail("integer literal out of range"); + tok.value = static_cast(v); + tok.isInt = true; + } + if (end == begin) + return fail("invalid numeric literal"); + const std::size_t len = end - begin; + // A hex literal whose final digit is e/E followed by '+' or '-' forms + // one single (invalid) pp-number in C++: cling rejects "0x1e+2" rather + // than computing 0x1e + 2. TFormula strips whitespace before compiling, + // so "0x1e + 2" is equally invalid. Keep such formulas failing. + if (isHex && (end[-1] == 'e' || end[-1] == 'E')) { + std::size_t j = i + len; + while (j < n && std::isspace(static_cast(_s[j]))) + ++j; + if (j < n && (_s[j] == '+' || _s[j] == '-')) + return fail("hex literal followed by an exponent-like sign is invalid in C++"); + } + // Reject trailing characters that would make this an invalid or + // differently-typed literal in C++ (e.g. "1e", "08", "1.5f", "1.2.3"). + if (i + len < n) { + const char next = _s[i + len]; + if (isIdentChar(next) || next == '.') + return fail("invalid numeric literal"); + } + tok.text = std::string_view{begin, len}; + out.push_back(tok); + i += len; + return true; + } + + /// Lex an identifier (including `::`-qualified names like TMath::Erf as one + /// token), or an `x[i]` variable reference. Only the exact shape + /// `x[]` is a variable; RooFormula::processFormula() guarantees it. + bool lexIdentOrVar(std::size_t &i, std::vector &out) + { + const std::size_t n = _s.size(); + const std::size_t start = i; + while (i < n && isIdentChar(_s[i])) + ++i; + // absorb `::name` qualifications (a single ':' is a ternary colon) + while (i + 2 < n && _s[i] == ':' && _s[i + 1] == ':' && isIdentStart(_s[i + 2])) { + i += 2; + while (i < n && isIdentChar(_s[i])) + ++i; + } + std::string_view name{_s.c_str() + start, i - start}; + + if (name == "x" && i < n && _s[i] == '[') { + std::size_t j = i + 1; + std::uint32_t index = 0; + if (j >= n || !std::isdigit(static_cast(_s[j]))) + return fail("malformed variable reference"); + while (j < n && std::isdigit(static_cast(_s[j]))) { + index = index * 10 + (_s[j] - '0'); + if (index > 1000000) + return fail("variable index out of range"); + ++j; + } + if (j >= n || _s[j] != ']') + return fail("malformed variable reference"); + Token tok; + tok.kind = Tok::Var; + tok.varIndex = index; + out.push_back(tok); + i = j + 1; + return true; + } + + Token tok; + tok.kind = Tok::Ident; + tok.text = name; + out.push_back(tok); + return true; + } + + /// Lex one operator or punctuation token, with maximal munch for the + /// two-character operators. Anything not in the supported set (notably `%`, + /// single `&`, `|`, `=`, brackets outside `x[i]`, string literals) fails. + bool lexOperator(std::size_t &i, std::vector &out) + { + const std::size_t n = _s.size(); + const char c = _s[i]; + const char c2 = i + 1 < n ? _s[i + 1] : '\0'; + Tok kind; + std::size_t len = 1; + switch (c) { + case '+': kind = Tok::Plus; break; + case '-': kind = Tok::Minus; break; + case '*': + if (c2 == '*') { // TFormula rewrites `**` to `^` + kind = Tok::Caret; + len = 2; + } else { + kind = Tok::Star; + } + break; + case '/': kind = Tok::Slash; break; + case '^': kind = Tok::Caret; break; + case '(': kind = Tok::LParen; break; + case ')': kind = Tok::RParen; break; + case ',': kind = Tok::Comma; break; + case '?': kind = Tok::Question; break; + case ':': kind = Tok::Colon; break; + case '<': + if (c2 == '=') { + kind = Tok::Le; + len = 2; + } else { + kind = Tok::Lt; + } + break; + case '>': + if (c2 == '=') { + kind = Tok::Ge; + len = 2; + } else { + kind = Tok::Gt; + } + break; + case '!': + if (c2 == '=') { + kind = Tok::Ne; + len = 2; + } else { + kind = Tok::Not; + } + break; + case '=': + if (c2 != '=') + return fail("unsupported operator '='"); + kind = Tok::Eq; + len = 2; + break; + case '&': + if (c2 != '&') + return fail("unsupported operator '&'"); + kind = Tok::AndAnd; + len = 2; + break; + case '|': + if (c2 != '|') + return fail("unsupported operator '|'"); + kind = Tok::OrOr; + len = 2; + break; + default: return fail(std::string{"unsupported character '"} + c + "'"); + } + Token tok; + tok.kind = kind; + out.push_back(tok); + i += len; + return true; + } + + std::string const &_s; + std::string &_error; +}; + +/// Binary operator precedence, mirroring C++ exactly (what cling compiled). +/// Loosest binds first; `?:` sits below level 1 (handled in parseTernary) and +/// unary +/-/! above level 6 (handled in parseUnary), with `^`/`**` +/// exponentiation tighter still (handled in parsePower). All these operators +/// are left-associative. Bitwise `& | ^`(C++ meaning) are not part of the +/// dialect: `^` is exponentiation and single `&`/`|` fail to tokenize. +struct BinOpInfo { + Tok tok; + int prec; + Op op; +}; + +// clang-format off +constexpr BinOpInfo gBinaryOps[] = { + {Tok::OrOr, 1, Op::Or}, + {Tok::AndAnd, 2, Op::And}, + {Tok::Eq, 3, Op::EQ}, + {Tok::Ne, 3, Op::NE}, + {Tok::Lt, 4, Op::LT}, + {Tok::Le, 4, Op::LE}, + {Tok::Gt, 4, Op::GT}, + {Tok::Ge, 4, Op::GE}, + {Tok::Plus, 5, Op::Add}, + {Tok::Minus, 5, Op::Sub}, + {Tok::Star, 6, Op::Mul}, + {Tok::Slash, 6, Op::Div}, +}; +// clang-format on + +BinOpInfo const *findBinOp(Tok kind) +{ + for (auto const &info : gBinaryOps) { + if (info.tok == kind) + return &info; + } + return nullptr; +} + +class Parser { +public: + Parser(std::vector const &tokens, unsigned int nVars, std::string &error) + : _tokens{tokens}, _nVars{nVars}, _error{error} + { + } + + /// C++ int-vs-double typing of a subexpression, tracked to reject + /// constructs whose cling semantics double arithmetic cannot reproduce. + struct ExprInfo { + bool isInt = false; + }; + + bool run(Program &prog) + { + _used.assign(_nVars, false); + ExprInfo info; + if (!parseTernary(info)) + return false; + if (peek().kind != Tok::End) + return fail("unexpected token after end of expression"); + + // Compute the maximum evaluation stack depth. The switch is exhaustive + // over ExprOp with no default case, so that adding an opcode without + // extending the accounting is a compiler warning (-Wswitch); the + // static_assert additionally breaks the build when the enum grows. + static_assert(static_cast(Op::Call4) == 23, + "ExprOp changed: update the stack-depth accounting switch in RooFormulaParser"); + int depth = 0; + int maxDepth = 0; + for (Instr const &ins : _code) { + switch (ins.op) { + case Op::Const: + case Op::Var: ++depth; break; + case Op::Neg: + case Op::Not: + case Op::Sq: + case Op::IntNorm: + case Op::Call1: break; + case Op::Add: + case Op::Sub: + case Op::Mul: + case Op::Div: + case Op::LT: + case Op::LE: + case Op::GT: + case Op::GE: + case Op::EQ: + case Op::NE: + case Op::And: + case Op::Or: + case Op::Pow: + case Op::Call2: --depth; break; + case Op::Select: + case Op::Call3: depth -= 2; break; + case Op::Call4: depth -= 3; break; + } + maxDepth = std::max(maxDepth, depth); + } + if (maxDepth > static_cast(RooExprEvaluator::kMaxStackDepth)) + return fail("expression too deep"); + + prog.code = std::move(_code); + prog.stackDepth = maxDepth; + // Trim to the highest used index so that programs can be shared between + // formulas with different (sufficiently long) variable lists. + std::size_t lastUsed = 0; + for (std::size_t i = 0; i < _used.size(); ++i) { + if (_used[i]) + lastUsed = i + 1; + } + _used.resize(lastUsed); + prog.usedVars = std::move(_used); + return true; + } + +private: + static constexpr int kMaxRecursionDepth = 128; + + struct DepthGuard { + DepthGuard(int &d) : _d{d} { ++_d; } + ~DepthGuard() { --_d; } + int &_d; + }; + + bool fail(std::string msg) + { + if (_error.empty()) + _error = std::move(msg); + return false; + } + + Token const &peek() const { return _tokens[_pos]; } + Token const &next() { return _tokens[_pos++]; } + + void emit(Op op, std::uint32_t arg = 0, double konst = 0.0) + { + Instr ins; + ins.op = op; + ins.arg = arg; + ins.konst = konst; + _code.push_back(ins); + } + + /// conditional-expression: right-associative, both branches always + /// evaluated with the value of the active branch selected (Op::Select). + bool parseTernary(ExprInfo &out) + { + DepthGuard guard{_depth}; + if (_depth > kMaxRecursionDepth) + return fail("expression too deeply nested"); + if (!parseBinary(1, out)) + return false; + if (peek().kind != Tok::Question) + return true; + next(); + ExprInfo left; + ExprInfo right; + if (!parseTernary(left)) + return false; + if (peek().kind != Tok::Colon) + return fail("expected ':' in conditional expression"); + next(); + if (!parseTernary(right)) + return false; + emit(Op::Select); + // The C++ type of `c ? a : b` is int only if both branches are int. + out.isInt = left.isInt && right.isInt; + return true; + } + + /// Precedence climbing over gBinaryOps. + bool parseBinary(int minPrec, ExprInfo &out) + { + if (!parseUnary(out)) + return false; + while (true) { + BinOpInfo const *info = findBinOp(peek().kind); + if (!info || info->prec < minPrec) + break; + next(); + ExprInfo rhs; + if (!parseBinary(info->prec + 1, rhs)) + return false; + switch (info->op) { + case Op::Add: + case Op::Sub: + case Op::Mul: + out.isInt = out.isInt && rhs.isInt; + emit(info->op); + // cling would compute `int * int` in integer arithmetic, where + // e.g. (-1) * 0 is +0 and not the -0.0 of double arithmetic. + if (out.isInt && info->op == Op::Mul) + emit(Op::IntNorm); + break; + case Op::Div: + // `1/2` is a truncating integer division in cling. Not supported; + // fall back to TFormula so the behavior is unchanged. + if (out.isInt && rhs.isInt) + return fail("integer division has truncating semantics in TFormula/cling"); + out.isInt = false; + emit(info->op); + break; + default: + // comparisons and logical operators: C++ result type is bool + out.isInt = true; + emit(info->op); + break; + } + } + return true; + } + + /// unary-expression: prefix `+`, `-`, `!`. + bool parseUnary(ExprInfo &out) + { + DepthGuard guard{_depth}; + if (_depth > kMaxRecursionDepth) + return fail("expression too deeply nested"); + switch (peek().kind) { + case Tok::Plus: next(); return parseUnary(out); // unary plus: no-op, type preserved + case Tok::Minus: + next(); + if (!parseUnary(out)) + return false; + emit(Op::Neg); + if (out.isInt) + emit(Op::IntNorm); // cling: -(int)0 is +0, not -0.0 + return true; + case Tok::Not: + next(); + if (!parseUnary(out)) + return false; + emit(Op::Not); + out.isInt = true; + return true; + default: return parsePower(out); + } + } + + /// Exponentiation via `^` (or `**`), reproducing TFormula's textual + /// `a^b` -> `pow(a,b)` rewrite (TFormula::HandleExponentiation): + /// right-associative, binding tighter than `*`, `/` and unary minus, with + /// TFormula's special case `expr^2` -> `TMath::Sq(expr)` when the exponent + /// is spelled exactly `2`. + bool parsePower(ExprInfo &out) + { + const std::size_t startTok = _pos; + if (!parsePrimary(out)) + return false; + if (peek().kind != Tok::Caret) + return true; + // TFormula's textual operand scan runs through `,` and `:` (it only + // stops at operators and parentheses), so a `^` directly adjacent to a + // function-argument or ternary boundary produces invalid code today + // (e.g. `pow(x,2^3)` or `c?a:b^2` do not compile). Keep those failing. + if (startTok > 0 && (_tokens[startTok - 1].kind == Tok::Comma || _tokens[startTok - 1].kind == Tok::Colon)) { + return fail("'^' operand adjacent to ',' or ':' is invalid in TFormula"); + } + next(); + const std::size_t codeSizeBeforeRhs = _code.size(); + ExprInfo rhs; + bool rhsIsLiteralTwo = false; + if (!parsePowerRhs(rhs, rhsIsLiteralTwo)) + return false; + if (peek().kind == Tok::Comma || peek().kind == Tok::Colon) + return fail("'^' operand adjacent to ',' or ':' is invalid in TFormula"); + if (rhsIsLiteralTwo) { + // TMath::Sq is unary: drop the emitted `Const 2` again. + _code.resize(codeSizeBeforeRhs); + emit(Op::Sq); + } else { + emit(Op::Pow); + } + out.isInt = false; // pow() and TMath::Sq(Double_t) return double + return true; + } + + /// The right-hand side of `^`: an optional single sign, then a + /// power-expression (making `^` right-associative, `x^-2^3` = pow(x,-(2^3))). + bool parsePowerRhs(ExprInfo &out, bool &isLiteralTwo) + { + DepthGuard guard{_depth}; + if (_depth > kMaxRecursionDepth) + return fail("expression too deeply nested"); + bool haveSign = false; + bool negate = false; + if (peek().kind == Tok::Plus) { + next(); + haveSign = true; + } else if (peek().kind == Tok::Minus) { + next(); + haveSign = true; + negate = true; + } + const std::size_t operandStart = _pos; + if (!parsePower(out)) + return false; + // TFormula turns `a^b` into TMath::Sq(a) only when the exponent is the + // literal token `2` (not `2.0`, not `(2)`, not `+2`). + isLiteralTwo = !haveSign && _pos == operandStart + 1 && _tokens[operandStart].kind == Tok::Number && + _tokens[operandStart].text == "2"; + if (negate) { + emit(Op::Neg); + if (out.isInt) + emit(Op::IntNorm); + } + return true; + } + + /// primary-expression: literal, `x[i]`, parenthesized expression, or an + /// allow-listed function call (zero-argument constant calls are folded). + bool parsePrimary(ExprInfo &out) + { + switch (peek().kind) { + case Tok::Number: { + Token const &tok = next(); + emit(Op::Const, 0, tok.value); + out.isInt = tok.isInt; + return true; + } + case Tok::Var: { + Token const &tok = next(); + if (tok.varIndex >= _nVars) + return fail("formula references x[" + std::to_string(tok.varIndex) + "] but fewer variables were provided"); + _used[tok.varIndex] = true; + emit(Op::Var, tok.varIndex); + out.isInt = false; + return true; + } + case Tok::LParen: { + next(); + if (!parseTernary(out)) + return false; + if (peek().kind != Tok::RParen) + return fail("expected ')'"); + next(); + return true; + } + case Tok::Ident: return parseCall(out); + default: return fail("expected an expression"); + } + } + + bool parseCall(ExprInfo &out) + { + const std::string name{next().text}; + if (peek().kind != Tok::LParen) + return fail("unknown identifier '" + name + "'"); + next(); + unsigned int nArgs = 0; + ExprInfo argInfo[4]; + if (peek().kind == Tok::RParen) { + next(); + } else { + while (true) { + if (nArgs == 4) + return fail("too many arguments in call to '" + name + "'"); + if (!parseTernary(argInfo[nArgs])) + return false; + ++nArgs; + if (peek().kind == Tok::Comma) { + next(); + continue; + } + if (peek().kind == Tok::RParen) { + next(); + break; + } + return fail("expected ',' or ')' in call to '" + name + "'"); + } + } + + std::uint32_t index = 0; + RooFormulaFunctions::Entry const *entry = RooFormulaFunctions::find(name, nArgs, index); + if (!entry) { + return fail("unsupported function '" + name + "' with " + std::to_string(nArgs) + " argument(s)"); + } + + using RooFormulaFunctions::TypeRule; + switch (entry->rule) { + case TypeRule::Double: out.isInt = false; break; + case TypeRule::SameAsFirstArg: out.isInt = argInfo[0].isInt; break; + case TypeRule::Int: out.isInt = true; break; + case TypeRule::MinMax: + // e.g. std::min(x, 3) with double x and int 3 does not compile in + // cling, so such formulas are invalid in TFormula today. Keep it so. + if (argInfo[0].isInt != argInfo[1].isInt) + return fail("'" + name + "' with mixed int/double arguments is invalid in TFormula"); + out.isInt = argInfo[0].isInt; + break; + } + + switch (nArgs) { + case 0: + // Zero-argument calls are the TMath constants: fold to a literal. + emit(Op::Const, 0, entry->fn0()); + break; + case 1: emit(Op::Call1, index); break; + case 2: emit(Op::Call2, index); break; + case 3: emit(Op::Call3, index); break; + case 4: emit(Op::Call4, index); break; + } + if (out.isInt) + emit(Op::IntNorm); // int-valued calls cannot yield -0.0 in cling + return true; + } + + std::vector const &_tokens; + unsigned int _nVars = 0; + std::string &_error; + std::size_t _pos = 0; + int _depth = 0; + std::vector _code; + std::vector _used; +}; + +std::shared_ptr parseImpl(std::string const &formula, unsigned int nVars, std::string &error) +{ + std::vector tokens; + if (!Tokenizer{formula, error}.run(tokens)) + return nullptr; + auto prog = std::make_shared(); + if (!Parser{tokens, nVars, error}.run(*prog)) + return nullptr; + prog->formula = formula; + return prog; +} + +} // namespace + +std::shared_ptr +RooFormulaParser::compile(std::string const &processedFormula, unsigned int nVars, std::string *error) +{ + // Process-wide registry of compiled programs, keyed on the processed + // formula string: identical formulas (e.g. thousands of structurally equal + // HistFactory expressions) share one immutable instruction vector. Like + // TFormula's gClingFunctions cache, the registry grows without bound over + // the process lifetime; entries are small (the instruction vector). + // Only successful parses are cached: formulas destined for the TFormula + // fallback are re-parsed on each construction, which is cheap compared to + // the JIT compilation that follows. + static std::mutex mutex; + static std::unordered_map> registry; + + std::string errorBuffer; + std::string &err = error ? *error : errorBuffer; + + { + std::lock_guard lock{mutex}; + auto it = registry.find(processedFormula); + if (it != registry.end()) { + if (it->second->usedVars.size() <= nVars) + return it->second; + err = "formula references more variables than provided"; + return nullptr; + } + } + + auto prog = parseImpl(processedFormula, nVars, err); + if (!prog) + return nullptr; + + std::lock_guard lock{mutex}; + // If another thread compiled the same formula concurrently, share its copy. + return registry.emplace(processedFormula, std::move(prog)).first->second; +} + +/// \endcond diff --git a/roofit/roofitcore/src/RooFormulaParser.h b/roofit/roofitcore/src/RooFormulaParser.h new file mode 100644 index 0000000000000..a8c228ce46a11 --- /dev/null +++ b/roofit/roofitcore/src/RooFormulaParser.h @@ -0,0 +1,44 @@ +/// \cond ROOFIT_INTERNAL + +/* + * Project: RooFit + * + * Copyright (c) 2026, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +#ifndef ROO_FORMULA_PARSER +#define ROO_FORMULA_PARSER + +#include "RooExprEvaluator.h" + +#include +#include + +namespace RooFormulaParser { + +/// Try to compile the given processed formula string (all variables in the +/// `x[i]` dialect produced by RooFormula::processFormula()) into a program for +/// the JIT-free RooExprEvaluator. +/// +/// Returns nullptr on *any* construct the JIT-free evaluator does not support +/// (unknown identifier, unknown function, unsupported operator, ...); the +/// caller then falls back to the TFormula backend. If `error` is non-null, it +/// is filled with the reason on failure. +/// +/// Successfully compiled programs are cached in a process-wide registry keyed +/// on the formula string, so identical formulas share one immutable +/// instruction vector (mirroring what TFormula's gClingFunctions cache does +/// for JIT-compiled code). Compilation takes a mutex; evaluation of the +/// returned program is lock-free. +std::shared_ptr +compile(std::string const &processedFormula, unsigned int nVars, std::string *error = nullptr); + +} // namespace RooFormulaParser + +#endif + +/// \endcond diff --git a/roofit/roofitcore/src/RooFormulaUtils.cxx b/roofit/roofitcore/src/RooFormulaUtils.cxx index 2ccb69b25d2ac..461910968aea7 100644 --- a/roofit/roofitcore/src/RooFormulaUtils.cxx +++ b/roofit/roofitcore/src/RooFormulaUtils.cxx @@ -24,6 +24,15 @@ Free functions to translate and evaluate user-defined expressions of RooAbsArgs. See RooFormulaUtils.h for a description of the supported expression dialect. To debug the formula preprocessing, activate the RooFit::DEBUG message level for the RooFit::InputArguments topic. + +### Evaluation backends +By default, makeEvaluator() compiles the expression with a small built-in +parser and evaluates it without any use of the interpreter/JIT. Expressions +the parser does not support silently fall back to the traditional TFormula +(cling JIT) backend, so any expression that worked before keeps working. The +environment variable `ROOFIT_FORMULA_BACKEND` overrides this: `tformula` +always uses the TFormula backend, and `ast` disables the fallback, turning +unsupported expressions into hard errors (useful for testing). **/ #include "RooFormulaUtils.h" @@ -35,14 +44,18 @@ RooFit::DEBUG message level for the RooFit::InputArguments topic. #include "RooCurve.h" #include "RooFitImplHelpers.h" #include "RooMsgService.h" +#include "RooExprEvaluator.h" +#include "RooFormulaParser.h" #include "RooTFormulaEvaluator.h" #include "TFormula.h" #include #include +#include #include #include +#include #include #include @@ -50,6 +63,65 @@ using std::sregex_iterator; namespace { +/// Which evaluation backend the formula should use, from ROOFIT_FORMULA_BACKEND. +enum class FormulaBackend { + AstWithFallback = 0, ///< default: JIT-free parser, silent TFormula fallback + AstOnly = 1, ///< `ast`: fail loudly instead of falling back + TFormulaOnly = 2 ///< `tformula`: always use the TFormula backend +}; + +std::mutex gFormulaBackendMutex; +int gFormulaBackendCached = -1; ///< -1: environment variable not read yet +bool gFormulaBackendWarnedChange = false; +std::string gFormulaBackendEnvSeen; ///< raw env value at first read ("" if unset) + +/// Read ROOFIT_FORMULA_BACKEND once (thread-safely) and cache the result. On +/// later reads, warn (once) if the environment variable no longer matches the +/// cached value: changing it after the first evaluator creation has no +/// effect. +FormulaBackend formulaBackend() +{ + std::lock_guard lock{gFormulaBackendMutex}; + const char *env = std::getenv("ROOFIT_FORMULA_BACKEND"); + const std::string val = env ? env : ""; + if (gFormulaBackendCached < 0) { + FormulaBackend mode = FormulaBackend::AstWithFallback; + if (val == "ast") { + mode = FormulaBackend::AstOnly; + } else if (val == "tformula") { + mode = FormulaBackend::TFormulaOnly; + } else if (!val.empty()) { + oocoutW(nullptr, InputArguments) << "Ignoring unknown ROOFIT_FORMULA_BACKEND value '" << val + << "' (supported: \"ast\", \"tformula\")" << std::endl; + } + gFormulaBackendCached = static_cast(mode); + gFormulaBackendEnvSeen = val; + gFormulaBackendWarnedChange = false; + } else if (val != gFormulaBackendEnvSeen && !gFormulaBackendWarnedChange) { + gFormulaBackendWarnedChange = true; + oocoutW(nullptr, InputArguments) << "ROOFIT_FORMULA_BACKEND changed from '" << gFormulaBackendEnvSeen << "' to '" + << val << "' after it was first read; the change has no effect in this process" + << std::endl; + } + return static_cast(gFormulaBackendCached); +} + +} // namespace + +namespace RooFormulaInternal { + +void resetFormulaBackendForTesting() +{ + std::lock_guard lock{gFormulaBackendMutex}; + gFormulaBackendCached = -1; + gFormulaBackendWarnedChange = false; + gFormulaBackendEnvSeen.clear(); +} + +} // namespace RooFormulaInternal + +namespace { + /// Convert `@i`-style references to `x[i]`. void convertArobaseReferences(std::string &formula) { @@ -328,6 +400,16 @@ RooFormulaUtils::reconstructFormula(std::string internalRepr, RooArgList const & /// Create the evaluation engine for a processed formula, checking that the /// formula compiles and also fulfills the assumptions. Throws on failure, /// with the original formula string appearing in the error messages. +/// +/// First, the JIT-free expression parser is tried (unless disabled via +/// ROOFIT_FORMULA_BACKEND=tformula). On any unsupported construct it silently +/// falls back to the TFormula backend, so genuinely invalid formulas produce +/// exactly the same errors as before. +/// \param[in] name Name of the calling object, used to name the engine and in error messages. +/// \param[in] processedFormula The formula string in the normalized `x[i]` dialect, +/// with `i` referring to the position in `varList`. +/// \param[in] origFormula The original formula string as given by the user, used in error messages. +/// \param[in] varList List of variables to be passed to the formula. std::unique_ptr RooFormulaUtils::makeEvaluator(std::string const &name, std::string const &processedFormula, std::string const &origFormula, RooArgList const &varList) @@ -337,6 +419,40 @@ RooFormulaUtils::makeEvaluator(std::string const &name, std::string const &proce << "\n\t" << processedFormula << "\n and used as" << "\n\t" << reconstructFormula(processedFormula, varList) << "\n with the parameters " << varList << std::endl; + const FormulaBackend backend = formulaBackend(); + + if (backend != FormulaBackend::TFormulaOnly) { + std::string parseError; + if (auto program = RooFormulaParser::compile(processedFormula, varList.size(), &parseError)) { + return std::make_unique(std::move(program)); + } + if (backend == FormulaBackend::AstOnly) { + std::stringstream msg; + msg << "RooFormula '" << name << "' could not be compiled by the RooFit expression parser (" << parseError + << "), and ROOFIT_FORMULA_BACKEND=ast disables the TFormula fallback." + << "\nInput:\n\t" << origFormula << "\nProcessed formula:\n\t" << processedFormula << std::endl; + oocoutF(static_cast(nullptr), InputArguments) << msg.str(); + throw std::runtime_error(msg.str()); + } + // Report the silent fallback once per process at INFO level (further + // fallbacks are only visible on the debug stream, to avoid spamming). + static std::once_flag fallbackNoticeFlag; + std::call_once(fallbackNoticeFlag, [&] { + oocoutI(static_cast(nullptr), InputArguments) + << "RooFormula '" << name + << "': expression not supported by the RooFit formula " + "parser (" + << parseError + << "), falling back to the TFormula (cling JIT) backend. This notice is only " + "printed once; see the ROOFIT_FORMULA_BACKEND environment variable and the InputArguments debug " + "stream for details." + << std::endl; + }); + oocxcoutD(static_cast(nullptr), InputArguments) + << "RooFormula '" << name << "': expression not supported by the RooFit expression parser (" << parseError + << "), falling back to TFormula" << std::endl; + } + return std::make_unique(name.c_str(), processedFormula, origFormula, varList); } diff --git a/roofit/roofitcore/src/RooFormulaUtils.h b/roofit/roofitcore/src/RooFormulaUtils.h index e792af4a950d0..6330bdc9b3889 100644 --- a/roofit/roofitcore/src/RooFormulaUtils.h +++ b/roofit/roofitcore/src/RooFormulaUtils.h @@ -83,6 +83,14 @@ std::list *plotSamplingHint(BinningMap const &binnings, RooArgList const } // namespace RooFormulaUtils +namespace RooFormulaInternal { + +/// Testing hook: discard the cached ROOFIT_FORMULA_BACKEND setting so that it +/// is read again from the environment on the next evaluator creation. +void resetFormulaBackendForTesting(); + +} // namespace RooFormulaInternal + #endif /// \endcond diff --git a/roofit/roofitcore/src/RooFormulaVar.cxx b/roofit/roofitcore/src/RooFormulaVar.cxx index ca740c48e6dc9..c19c1c2e019c4 100644 --- a/roofit/roofitcore/src/RooFormulaVar.cxx +++ b/roofit/roofitcore/src/RooFormulaVar.cxx @@ -326,9 +326,30 @@ double RooFormulaVar::defaultErrorLevel() const return 1.0 ; } +//////////////////////////////////////////////////////////////////////////////// +/// Name of the cling-JIT-compiled function that evaluates this formula, which +/// generated code from the codegen path calls by name. +/// +/// If the formula is evaluated by the TFormula backend, this is the function +/// of the evaluating TFormula itself. Otherwise, a TFormula is created lazily, +/// only to serve the codegen path. +/// +/// TODO(Phase 2.5): remove the lazily-created TFormula once codegen emits C++ +/// for the expression directly from the parsed representation instead of +/// calling the JIT-compiled TFormula function. Until then, codegen keeps +/// working exactly as before, at the cost of one JIT compilation per formula +/// -- but only when codegen is actually used. Like codegen itself, this lazy +/// creation is not thread-safe. std::string RooFormulaVar::getUniqueFuncName() const { - return evaluator().getTFormula()->GetUniqueFuncName().Data(); + if (TFormula *tFormula = evaluator().getTFormula()) { + return tFormula->GetUniqueFuncName().Data(); + } + if (!_tFormulaForCodegen) { + // evaluator() above has normalized _formExpr to the processed `x[i]` dialect. + _tFormulaForCodegen = std::make_unique(GetName(), _formExpr.Data(), /*addToGlobList=*/false); + } + return _tFormulaForCodegen->GetUniqueFuncName().Data(); } std::unique_ptr diff --git a/roofit/roofitcore/src/RooGenericPdf.cxx b/roofit/roofitcore/src/RooGenericPdf.cxx index 549c293995f1b..ea411c6706a3f 100644 --- a/roofit/roofitcore/src/RooGenericPdf.cxx +++ b/roofit/roofitcore/src/RooGenericPdf.cxx @@ -227,7 +227,20 @@ void RooGenericPdf::writeToStream(ostream& os, bool compact) const } } +//////////////////////////////////////////////////////////////////////////////// +/// Name of the cling-JIT-compiled function that evaluates this formula, which +/// generated code from the codegen path calls by name. +/// +/// See RooFormulaVar::getUniqueFuncName() for the details of the lazily +/// created TFormula on the JIT-free expression backend. std::string RooGenericPdf::getUniqueFuncName() const { - return evaluator().getTFormula()->GetUniqueFuncName().Data(); + if (TFormula *tFormula = evaluator().getTFormula()) { + return tFormula->GetUniqueFuncName().Data(); + } + if (!_tFormulaForCodegen) { + // evaluator() above has normalized _formExpr to the processed `x[i]` dialect. + _tFormulaForCodegen = std::make_unique(GetName(), _formExpr.Data(), /*addToGlobList=*/false); + } + return _tFormulaForCodegen->GetUniqueFuncName().Data(); } diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index ac9d5bb73b7f0..b0c0d0da706d7 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -64,6 +64,7 @@ ROOT_ADD_GTEST(testRooAbsPdf testRooAbsPdf.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooDataSet testRooDataSet.cxx LIBRARIES Tree RooFitCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/dataSet_with_errors_6_26_10.root) ROOT_ADD_GTEST(testRooFormula testRooFormula.cxx LIBRARIES RooFitCore ROOT::TestSupport) +ROOT_ADD_GTEST(testRooFormulaEvaluator testRooFormulaEvaluator.cxx LIBRARIES RooFitCore Hist MathCore) ROOT_ADD_GTEST(testProxiesAndCategories testProxiesAndCategories.cxx LIBRARIES RooFitCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testProxiesAndCategories_1.root diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx new file mode 100644 index 0000000000000..0b4cbc27e595b --- /dev/null +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -0,0 +1,848 @@ +// Tests for the JIT-free RooFit formula evaluation backend +// (RooFormulaParser + RooExprEvaluator), and its silent-TFormula-fallback +// contract. +// Author: Jonas Rembser, CERN 2026 + +#include "../src/RooFormulaUtils.h" +#include "../src/RooFormulaParser.h" +#include "../src/RooExprEvaluator.h" + +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +/// Bitwise comparison; any two NaNs count as equal. +bool sameBits(double a, double b) +{ + if (std::isnan(a) && std::isnan(b)) + return true; + return std::memcmp(&a, &b, sizeof(double)) == 0; +} + +/// Compile with the JIT-free parser and evaluate. The expression must parse. +double astVal(std::string const &expr, std::vector const &vars = {}) +{ + auto prog = RooFormulaParser::compile(expr, vars.size()); + if (!prog) { + ADD_FAILURE() << "expression unexpectedly failed to parse: " << expr; + return std::numeric_limits::quiet_NaN(); + } + return RooExprEvaluator{prog}.eval(vars.data()); +} + +bool astParses(std::string const &expr, unsigned int nVars = 1) +{ + return RooFormulaParser::compile(expr, nVars) != nullptr; +} + +/// Whether the evaluation engine is the JIT-free expression backend. +bool isAstBackend(RooFormulaEvaluator const &ev) +{ + return dynamic_cast(&ev) != nullptr; +} + +/// Set or unset ROOFIT_FORMULA_BACKEND for the lifetime of the object, +/// resetting the backend's read-once cache on both ends. +class ScopedBackendEnv { +public: + ScopedBackendEnv(const char *value) + { + if (const char *old = std::getenv("ROOFIT_FORMULA_BACKEND")) + _old = old; + if (value) + setenv("ROOFIT_FORMULA_BACKEND", value, /*overwrite=*/1); + else + unsetenv("ROOFIT_FORMULA_BACKEND"); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + ~ScopedBackendEnv() + { + if (_old.empty()) + unsetenv("ROOFIT_FORMULA_BACKEND"); + else + setenv("ROOFIT_FORMULA_BACKEND", _old.c_str(), 1); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + +private: + std::string _old; +}; + +} // namespace + +// Precedence must match C++ exactly (that is what cling compiled). Each case +// evaluates to a different number under a plausible wrong precedence. +TEST(RooFormulaEvaluator, Precedence) +{ + EXPECT_EQ(astVal("1+2*3"), 7.0); + EXPECT_EQ(astVal("6-8/2."), 2.0); + EXPECT_EQ(astVal("2*3+4*5"), 26.0); + EXPECT_EQ(astVal("1-2-3"), -4.0); // left-associative + EXPECT_EQ(astVal("10-2+3"), 11.0); // left-associative + EXPECT_EQ(astVal("16/4./2."), 2.0); // left-associative + EXPECT_EQ(astVal("1<2==1"), 1.0); // (1<2)==1, not 1<(2==1) + EXPECT_EQ(astVal("1||0&&0"), 1.0); // 1||(0&&0), not (1||0)&&0 + EXPECT_EQ(astVal("1+1<1+3"), 1.0); // (1+1)<(1+3) + EXPECT_EQ(astVal("0<1==0<1"), 1.0); // relational binds tighter than equality + EXPECT_EQ(astVal("1?2:3+10"), 2.0); // ternary is loosest: 1?2:(13) + EXPECT_EQ(astVal("1?2:3?4:5"), 2.0); // right-associative + EXPECT_EQ(astVal("0?2:3?4:5"), 4.0); // 0?2:(3?4:5) + EXPECT_EQ(astVal("0?2:0?4:5"), 5.0); + EXPECT_EQ(astVal("1?0?6:7:8"), 7.0); // nested in the middle operand + EXPECT_EQ(astVal("2--3"), 5.0); // binary minus then unary minus + EXPECT_EQ(astVal("2- -3"), 5.0); + EXPECT_EQ(astVal("-2*3"), -6.0); + EXPECT_EQ(astVal("!0"), 1.0); + EXPECT_EQ(astVal("!3"), 0.0); + EXPECT_EQ(astVal("!!3"), 1.0); + EXPECT_EQ(astVal("!1<2"), 1.0); // (!1)<2 +} + +// `^` (and `**`) follow TFormula::HandleExponentiation: right-associative, +// tighter than `*`, `/` and unary minus, one optional sign in the exponent, +// and the textual exponent `2` becomes TMath::Sq. Reference values were +// checked against TFormula/cling directly. +TEST(RooFormulaEvaluator, PowerOperator) +{ + EXPECT_EQ(astVal("2^3"), 8.0); + EXPECT_EQ(astVal("2**3"), 8.0); + EXPECT_EQ(astVal("-2^2"), -4.0); // -(2^2), not (-2)^2 + EXPECT_EQ(astVal("2^3^2"), 512.0); // right-assoc: 2^(3^2) + EXPECT_EQ(astVal("5*2^3"), 40.0); // not (5*2)^3 + EXPECT_EQ(astVal("2^2*3"), 12.0); // not 2^(2*3) + EXPECT_EQ(astVal("2^3+1"), 9.0); + EXPECT_EQ(astVal("x[0]^-2", {2.0}), 0.25); // sign in the exponent + EXPECT_TRUE(sameBits(astVal("x[0]^-2.5", {2.0}), std::pow(2.0, -2.5))); + EXPECT_TRUE(sameBits(astVal("x[0]^+2.5", {2.0}), std::pow(2.0, 2.5))); + EXPECT_TRUE(sameBits(astVal("x[0]^-2^2", {2.0}), std::pow(2.0, -4.0))); // pow(x,-(2^2)) + EXPECT_TRUE(sameBits(astVal("x[0]^2", {1.3}), 1.3 * 1.3)); // TMath::Sq special case + EXPECT_TRUE(sameBits(astVal("x[0]**2", {1.3}), 1.3 * 1.3)); + EXPECT_TRUE(sameBits(astVal("x[0]^2.0", {1.3}), std::pow(1.3, 2.0))); + EXPECT_TRUE(sameBits(astVal("sin(x[0])^2", {0.7}), std::sin(0.7) * std::sin(0.7))); + EXPECT_EQ(astVal("0x64^2"), 10000.0); + EXPECT_TRUE(sameBits(astVal("x[0]^sin(x[0])", {2.0}), std::pow(2.0, std::sin(2.0)))); + EXPECT_EQ(astVal("2^x[0]*3", {3.0}), 24.0); +} + +TEST(RooFormulaEvaluator, ComparisonsYieldExactly0Or1) +{ + EXPECT_TRUE(sameBits(astVal("x[0]>2", {3.0}), 1.0)); + EXPECT_TRUE(sameBits(astVal("x[0]>2", {1.0}), 0.0)); + EXPECT_TRUE(sameBits(astVal("(x[0]>2)*3", {3.0}), 3.0)); + EXPECT_TRUE(sameBits(astVal("x[0]<=0.5", {0.5}), 1.0)); + EXPECT_TRUE(sameBits(astVal("x[0]==0.25", {0.25}), 1.0)); + EXPECT_TRUE(sameBits(astVal("x[0]!=0.25", {0.25}), 0.0)); + // comparisons with NaN are false, like in C++ + const double nan = std::numeric_limits::quiet_NaN(); + EXPECT_TRUE(sameBits(astVal("x[0]<1", {nan}), 0.0)); + EXPECT_TRUE(sameBits(astVal("x[0]>=1", {nan}), 0.0)); + EXPECT_TRUE(sameBits(astVal("x[0]==x[0]", {nan}), 0.0)); + // && and || convert their operands like C++ bool conversion (NaN is true) + EXPECT_TRUE(sameBits(astVal("x[0]&&1", {nan}), 1.0)); + EXPECT_TRUE(sameBits(astVal("0.5&&0"), 0.0)); + EXPECT_TRUE(sameBits(astVal("0.0||0.25"), 1.0)); +} + +TEST(RooFormulaEvaluator, TernarySelectsExactValue) +{ + // Both branches are evaluated (no branching, matching a future vectorized + // path), but the selected value is exact. + EXPECT_TRUE(sameBits(astVal("x[0]>0 ? log(x[0]) : -1", {-2.0}), -1.0)); + EXPECT_TRUE(sameBits(astVal("x[0]>0 ? log(x[0]) : -1", {2.0}), std::log(2.0))); + const double nan = std::numeric_limits::quiet_NaN(); + EXPECT_TRUE(sameBits(astVal("x[0] ? 1 : 2", {nan}), 1.0)); // NaN converts to true +} + +// Every allow-listed function in every accepted spelling, against the exact +// call the JIT-compiled code would have made. +TEST(RooFormulaEvaluator, Functions) +{ + const double v = 0.7311; + auto check1 = [&](const char *expr, double expected) { + EXPECT_TRUE(sameBits(astVal(std::string(expr) + "(x[0])", {v}), expected)) << expr; + }; + check1("sqrt", std::sqrt(v)); + check1("std::sqrt", std::sqrt(v)); + check1("TMath::Sqrt", TMath::Sqrt(v)); + check1("exp", std::exp(v)); + check1("std::exp", std::exp(v)); + check1("TMath::Exp", TMath::Exp(v)); + check1("log", std::log(v)); + check1("std::log", std::log(v)); + check1("TMath::Log", TMath::Log(v)); + check1("log10", std::log10(v)); + check1("std::log10", std::log10(v)); + check1("TMath::Log10", TMath::Log10(v)); + check1("sin", std::sin(v)); + check1("std::sin", std::sin(v)); + check1("TMath::Sin", TMath::Sin(v)); + check1("cos", std::cos(v)); + check1("std::cos", std::cos(v)); + check1("TMath::Cos", TMath::Cos(v)); + check1("tan", std::tan(v)); + check1("std::tan", std::tan(v)); + check1("TMath::Tan", TMath::Tan(v)); + check1("asin", std::asin(v)); + check1("std::asin", std::asin(v)); + check1("TMath::ASin", TMath::ASin(v)); + check1("acos", std::acos(v)); + check1("std::acos", std::acos(v)); + check1("TMath::ACos", TMath::ACos(v)); + check1("atan", std::atan(v)); + check1("std::atan", std::atan(v)); + check1("TMath::ATan", TMath::ATan(v)); + check1("sinh", std::sinh(v)); + check1("std::sinh", std::sinh(v)); + check1("TMath::SinH", TMath::SinH(v)); + check1("cosh", std::cosh(v)); + check1("std::cosh", std::cosh(v)); + check1("TMath::CosH", TMath::CosH(v)); + check1("tanh", std::tanh(v)); + check1("std::tanh", std::tanh(v)); + check1("TMath::TanH", TMath::TanH(v)); + check1("asinh", std::asinh(v)); + check1("std::asinh", std::asinh(v)); + check1("TMath::ASinH", TMath::ASinH(v)); + check1("atanh", std::atanh(v)); + check1("std::atanh", std::atanh(v)); + check1("TMath::ATanH", TMath::ATanH(v)); + check1("floor", std::floor(v)); + check1("std::floor", std::floor(v)); + check1("TMath::Floor", TMath::Floor(v)); + check1("ceil", std::ceil(v)); + check1("std::ceil", std::ceil(v)); + check1("TMath::Ceil", TMath::Ceil(v)); + check1("erf", std::erf(v)); + check1("std::erf", std::erf(v)); + check1("TMath::Erf", TMath::Erf(v)); + check1("erfc", std::erfc(v)); + check1("std::erfc", std::erfc(v)); + check1("TMath::Erfc", TMath::Erfc(v)); + check1("tgamma", std::tgamma(v)); + check1("std::tgamma", std::tgamma(v)); + check1("lgamma", std::lgamma(v)); + check1("std::lgamma", std::lgamma(v)); + check1("abs", std::fabs(v)); + check1("std::abs", std::fabs(v)); + check1("fabs", std::fabs(v)); + check1("std::fabs", std::fabs(v)); + check1("TMath::Abs", TMath::Abs(v)); + check1("sq", v * v); + check1("TMath::Sq", TMath::Sq(v)); + + // acosh needs an argument >= 1 + const double w = 1.5; + EXPECT_TRUE(sameBits(astVal("acosh(x[0])", {w}), std::acosh(w))); + EXPECT_TRUE(sameBits(astVal("std::acosh(x[0])", {w}), std::acosh(w))); + EXPECT_TRUE(sameBits(astVal("TMath::ACosH(x[0])", {w}), TMath::ACosH(w))); + + // int() is a C++ functional cast: truncation towards zero + EXPECT_EQ(astVal("int(2.7)"), 2.0); + EXPECT_EQ(astVal("int(-2.7)"), -2.0); + + // TMath::SignBit uses std::signbit (note: true for -0.0) + EXPECT_EQ(astVal("TMath::SignBit(x[0])", {-2.0}), 1.0); + EXPECT_EQ(astVal("TMath::SignBit(x[0])", {2.0}), 0.0); + EXPECT_EQ(astVal("TMath::SignBit(x[0])", {-0.0}), 1.0); + + // two-argument functions + EXPECT_TRUE(sameBits(astVal("pow(x[0],x[1])", {1.7, 2.5}), std::pow(1.7, 2.5))); + EXPECT_TRUE(sameBits(astVal("std::pow(x[0],x[1])", {1.7, 2.5}), std::pow(1.7, 2.5))); + EXPECT_TRUE(sameBits(astVal("TMath::Power(x[0],x[1])", {1.7, 2.5}), TMath::Power(1.7, 2.5))); + EXPECT_TRUE(sameBits(astVal("atan2(x[0],x[1])", {1.0, 2.0}), std::atan2(1.0, 2.0))); + EXPECT_TRUE(sameBits(astVal("std::atan2(x[0],x[1])", {1.0, 2.0}), std::atan2(1.0, 2.0))); + EXPECT_TRUE(sameBits(astVal("TMath::ATan2(x[0],x[1])", {1.0, 2.0}), TMath::ATan2(1.0, 2.0))); + // TMath::ATan2 differs from std::atan2 for x == -0.0; keep each spelling exact + EXPECT_TRUE(sameBits(astVal("TMath::ATan2(x[0],x[1])", {0.0, -0.0}), TMath::ATan2(0.0, -0.0))); + EXPECT_TRUE(sameBits(astVal("atan2(x[0],x[1])", {0.0, -0.0}), std::atan2(0.0, -0.0))); + EXPECT_TRUE(sameBits(astVal("fmod(x[0],x[1])", {7.5, 2.0}), std::fmod(7.5, 2.0))); + EXPECT_TRUE(sameBits(astVal("std::fmod(x[0],x[1])", {7.5, 2.0}), std::fmod(7.5, 2.0))); + EXPECT_TRUE(sameBits(astVal("min(x[0],x[1])", {1.0, 2.0}), 1.0)); + EXPECT_TRUE(sameBits(astVal("max(x[0],x[1])", {1.0, 2.0}), 2.0)); + EXPECT_TRUE(sameBits(astVal("min(2,3)"), 2.0)); + EXPECT_TRUE(sameBits(astVal("std::max(x[0],3.0)", {5.0}), 5.0)); + // std::min/max and TMath::Min/Max have opposite NaN behavior; both must be exact + const double nan = std::numeric_limits::quiet_NaN(); + EXPECT_TRUE(std::isnan(astVal("min(x[0],1.0)", {nan}))); // std::min(NaN, 1) = NaN + EXPECT_TRUE(sameBits(astVal("min(1.0,x[0])", {nan}), 1.0)); // std::min(1, NaN) = 1 + EXPECT_TRUE(sameBits(astVal("TMath::Min(x[0],1.0)", {nan}), 1.0)); + EXPECT_TRUE(std::isnan(astVal("TMath::Min(1.0,x[0])", {nan}))); + // sign resolves to TMath::Sign, which is std::copysign for doubles + EXPECT_TRUE(sameBits(astVal("sign(1.5,x[0])", {-3.0}), -1.5)); + EXPECT_TRUE(sameBits(astVal("sign(1.5,x[0])", {3.0}), 1.5)); + EXPECT_TRUE(sameBits(astVal("sign(1.5,x[0])", {-0.0}), -1.5)); // copysign semantics + EXPECT_TRUE(sameBits(astVal("TMath::Sign(1.5,x[0])", {-0.0}), -1.5)); + EXPECT_TRUE(sameBits(astVal("x[0]*sign(1.,x[0]+2.)", {-3.0}), 3.0)); + + // zero-argument constants + EXPECT_TRUE(sameBits(astVal("TMath::Pi()"), TMath::Pi())); + EXPECT_TRUE(sameBits(astVal("TMath::TwoPi()"), TMath::TwoPi())); + EXPECT_TRUE(sameBits(astVal("TMath::PiOver2()"), TMath::PiOver2())); + EXPECT_TRUE(sameBits(astVal("TMath::E()"), TMath::E())); + + // TMath::Gaus with 1 to 4 arguments (default args mean=0, sigma=1, norm=false) + EXPECT_TRUE(sameBits(astVal("TMath::Gaus(x[0])", {1.0}), TMath::Gaus(1.0))); + EXPECT_TRUE(sameBits(astVal("TMath::Gaus(x[0],2)", {1.0}), TMath::Gaus(1.0, 2))); + EXPECT_TRUE(sameBits(astVal("TMath::Gaus(x[0],2,3)", {1.0}), TMath::Gaus(1.0, 2, 3))); + EXPECT_TRUE(sameBits(astVal("TMath::Gaus(x[0],2,3,1)", {1.0}), TMath::Gaus(1.0, 2, 3, true))); +} + +TEST(RooFormulaEvaluator, Literals) +{ + EXPECT_EQ(astVal("1e2"), 100.0); + EXPECT_EQ(astVal("2e+3"), 2000.0); + EXPECT_TRUE(sameBits(astVal("0.2e-6"), 0.2e-6)); + EXPECT_TRUE(sameBits(astVal("-7.94004e+06"), -7.94004e+06)); + EXPECT_EQ(astVal(".5"), 0.5); + EXPECT_EQ(astVal("1."), 1.0); + EXPECT_TRUE(sameBits(astVal("3.360779"), 3.360779)); + EXPECT_EQ(astVal("0x64"), 100.0); // hex, like cling + EXPECT_EQ(astVal("010"), 8.0); // octal, like cling + EXPECT_EQ(astVal("0"), 0.0); +} + +TEST(RooFormulaEvaluator, Variables) +{ + EXPECT_EQ(astVal("x[0]", {3.0}), 3.0); + EXPECT_EQ(astVal("x[1]+2*x[0]", {3.0, 4.0}), 10.0); + EXPECT_EQ(astVal("x[2]", {0.0, 0.0, 7.0}), 7.0); + EXPECT_EQ(astVal("-x[0]", {3.0}), -3.0); + + // used-variable tracking + auto prog = RooFormulaParser::compile("x[1]+1", 3); + ASSERT_TRUE(prog); + RooExprEvaluator ev{prog}; + EXPECT_FALSE(ev.usesVariable(0)); + EXPECT_TRUE(ev.usesVariable(1)); + EXPECT_FALSE(ev.usesVariable(2)); + EXPECT_EQ(ev.processedFormula(), "x[1]+1"); +} + +// Everything the JIT-free path does not support must return nullptr from the +// parser (and thus silently fall back to the TFormula backend). +TEST(RooFormulaEvaluator, FallbackTriggers) +{ + // unknown identifiers (undefined variables land here too) + EXPECT_FALSE(astParses("y+1")); + EXPECT_FALSE(astParses("unknownFunc(x[0])")); + EXPECT_FALSE(astParses("pi")); // TFormula constants are not part of the dialect + EXPECT_FALSE(astParses("TMath::Pi")); // constants require the call syntax + EXPECT_FALSE(astParses("ROOT::Math::normal_pdf(x[0],1.,2.)")); + // wrong arity + EXPECT_FALSE(astParses("erf(x[0],1.0)")); + EXPECT_FALSE(astParses("TMath::Gaus(x[0],1,2,3,4)")); + // `%` does not compile on doubles in cling, so TFormula rejects it today + EXPECT_FALSE(astParses("x[0]%2")); + // integer division truncates in cling; not reproduced in double arithmetic + EXPECT_FALSE(astParses("1/2")); + EXPECT_FALSE(astParses("7/2+x[0]")); + EXPECT_FALSE(astParses("(x[0]>1)/2")); + EXPECT_FALSE(astParses("int(x[0])/2")); + EXPECT_FALSE(astParses("-1/2")); + // ... but other all-integer arithmetic is value-identical in doubles + EXPECT_EQ(astVal("7*3+2"), 23.0); + EXPECT_EQ(astVal("7./2"), 3.5); + EXPECT_EQ(astVal("7/2."), 3.5); + // min/max with mixed int/double arguments does not compile in cling + EXPECT_FALSE(astParses("min(x[0],3)")); + EXPECT_FALSE(astParses("max(3,x[0])")); + EXPECT_FALSE(astParses("TMath::Min(x[0],3)")); + EXPECT_TRUE(astParses("min(x[0],3.0)")); + // TFormula's textual `^` rewrite breaks next to `,` or `:`; both are + // invalid in TFormula today and must stay invalid (fall back) + EXPECT_FALSE(astParses("pow(x[0],2^3)")); + EXPECT_FALSE(astParses("x[0]>0?1:x[0]^2")); + EXPECT_TRUE(astParses("x[0]>0?1:(x[0]^2)")); + // TFormula parameters and parametrized shortcuts + EXPECT_FALSE(astParses("[0]+x[0]")); + EXPECT_FALSE(astParses("x[0]+[0]*0xaf")); + EXPECT_FALSE(astParses("gaus")); + EXPECT_FALSE(astParses("pol1")); + EXPECT_FALSE(astParses("gaus(0)+pol1(3)")); + // string literals, malformed numbers, stray characters + EXPECT_FALSE(astParses("\"abc\"")); + EXPECT_FALSE(astParses("1e")); + EXPECT_FALSE(astParses("1.5f")); + EXPECT_FALSE(astParses("08")); + EXPECT_FALSE(astParses("x[0] & 1")); + EXPECT_FALSE(astParses("x[0] | 1")); + EXPECT_FALSE(astParses("x[0] = 1")); + EXPECT_FALSE(astParses("")); + // malformed variable references + EXPECT_FALSE(astParses("x[a]")); + EXPECT_FALSE(astParses("x[0")); + EXPECT_FALSE(astParses("x[0+1]")); + // referencing more variables than provided + EXPECT_FALSE(RooFormulaParser::compile("x[1]+1", 1)); + EXPECT_TRUE(RooFormulaParser::compile("x[1]+1", 2)); + // deep nesting falls back instead of overflowing the parser stack + std::string deep(300, '('); + deep += "1"; + deep += std::string(300, ')'); + EXPECT_FALSE(astParses(deep)); + std::string shallow(50, '('); + shallow += "1"; + shallow += std::string(50, ')'); + EXPECT_TRUE(astParses(shallow)); +} + +// "0x1e+2" is one single (invalid) pp-number in C++, not 0x1e + 2: cling +// refused it, so formula construction threw. The lexer must not split +// it into two tokens; such formulas keep failing via the TFormula fallback. +// TFormula strips whitespace before compiling, so "0x1e + 2" is equally +// invalid. +TEST(RooFormulaEvaluator, HexLiteralWithExponentSign) +{ + EXPECT_FALSE(astParses("0x1e+2")); + EXPECT_FALSE(astParses("0x1E-2")); + EXPECT_FALSE(astParses("0x1e + 2")); + EXPECT_FALSE(astParses("x[0]+0x1e+2")); + EXPECT_FALSE(astParses("0x1e+x[0]")); + EXPECT_EQ(astVal("0x1e"), 30.0); + EXPECT_EQ(astVal("0x1f+2"), 33.0); // final digit not e/E: ordinary addition + EXPECT_EQ(astVal("2+0x1e"), 32.0); // sign before the literal is fine + EXPECT_EQ(astVal("x[0]-0x1e", {5.0}), -25.0); + EXPECT_EQ(astVal("(0x1e)+2"), 32.0); // ')' ends the pp-number + + // On the default backend the formula must still throw at construction, + // exactly like before (via the TFormula fallback path). + ScopedBackendEnv env{nullptr}; + RooRealVar x("x", "x", 1.0); + { + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TFormula::InputFormulaIntoCling", "Error compiling formula expression in Cling", + false); + diags.requiredDiag(kError, "TFormula::ProcessFormula", " is invalid", false); + diags.optionalDiag(kError, "prepareMethod", "Can't compile function TFormula", false); + diags.optionalDiag(kError, "cling", "invalid suffix", false); + EXPECT_THROW(RooFormulaUtils::makeFormulaEvaluator("f", "0x1e+2+x", RooArgList{x}), std::runtime_error); + } +} + +// A long chain of `^` recursed once per operator without a depth guard and +// segfaulted on parser stack overflow. It must fail the parse cleanly (and +// fall back) like deeply nested parentheses do. +TEST(RooFormulaEvaluator, DeepPowerChainFallsBack) +{ + std::string deep = "1"; + for (int i = 0; i < 200000; ++i) { + deep += "^1"; + } + EXPECT_FALSE(astParses(deep)); + // a modest chain still parses + std::string shallow = "2"; + for (int i = 0; i < 50; ++i) { + shallow += "^1"; + } + EXPECT_TRUE(astParses(shallow)); + EXPECT_EQ(astVal(shallow), 2.0); +} + +// Identical formula strings share one immutable program via the registry. +TEST(RooFormulaEvaluator, ProgramSharing) +{ + auto p1 = RooFormulaParser::compile("x[0]*2+sin(x[0])", 1); + auto p2 = RooFormulaParser::compile("x[0]*2+sin(x[0])", 5); + ASSERT_TRUE(p1); + EXPECT_EQ(p1.get(), p2.get()); +} + +TEST(RooFormulaEvaluator, RooFormulaIntegration) +{ + // This test is about the default backend; shield it from an ambient + // ROOFIT_FORMULA_BACKEND setting. + ScopedBackendEnv env{nullptr}; + + RooRealVar x("x", "x", 2.0); + RooRealVar y("y", "y", 3.0); + RooArgList vars{x, y}; + + auto f = RooFormulaUtils::makeFormulaEvaluator("f", "x*y+sin(x)", vars); + EXPECT_TRUE(isAstBackend(*f)); + EXPECT_TRUE(sameBits(RooFormulaUtils::evalFormula(*f, vars), 2.0 * 3.0 + std::sin(2.0))); + + // copies share the backend and evaluate identically + auto fCopy = f->clone(); + EXPECT_TRUE(isAstBackend(*fCopy)); + EXPECT_TRUE(sameBits(RooFormulaUtils::evalFormula(*fCopy, vars), RooFormulaUtils::evalFormula(*f, vars))); + + // used-variable pruning in the owning classes works on the AST path + RooFormulaVar g("g", "y*2", {x, y}); + EXPECT_EQ(g.dependents().size(), 1u); + EXPECT_TRUE(g.dependents().find("y")); + + // the stored expression is the processed formula (persistence!) + RooFormulaVar fVar("fVar", "x*y+sin(x)", {x, y}); + EXPECT_STREQ(fVar.expression(), "x[0]*x[1]+sin(x[0])"); + EXPECT_TRUE(sameBits(fVar.getVal(), 2.0 * 3.0 + std::sin(2.0))); + + // the codegen accessor lazily provides a JIT-compiled TFormula function + // even on the AST path + EXPECT_FALSE(fVar.getUniqueFuncName().empty()); +} + +TEST(RooFormulaEvaluator, BackendOverride) +{ + RooRealVar x("x", "x", 2.0); + + RooArgList vars{x}; + { + ScopedBackendEnv env{nullptr}; // default: AST with silent fallback + auto f = RooFormulaUtils::makeFormulaEvaluator("f", "x*2", vars); + EXPECT_TRUE(isAstBackend(*f)); + EXPECT_EQ(RooFormulaUtils::evalFormula(*f, vars), 4.0); + // unsupported expression: silently falls back and still works + auto g = RooFormulaUtils::makeFormulaEvaluator("g", "ROOT::Math::normal_pdf(x,1.,0.)", vars); + EXPECT_FALSE(isAstBackend(*g)); + EXPECT_TRUE(sameBits(RooFormulaUtils::evalFormula(*g, vars), ROOT::Math::normal_pdf(2.0, 1.0, 0.0))); + } + { + ScopedBackendEnv env{"tformula"}; + auto f = RooFormulaUtils::makeFormulaEvaluator("f", "x*2", vars); + EXPECT_FALSE(isAstBackend(*f)); + EXPECT_EQ(RooFormulaUtils::evalFormula(*f, vars), 4.0); + } + { + ScopedBackendEnv env{"ast"}; + auto f = RooFormulaUtils::makeFormulaEvaluator("f", "x*2", vars); + EXPECT_TRUE(isAstBackend(*f)); + EXPECT_EQ(RooFormulaUtils::evalFormula(*f, vars), 4.0); + // unsupported expression: fail loudly instead of falling back + EXPECT_THROW(RooFormulaUtils::makeFormulaEvaluator("g", "ROOT::Math::normal_pdf(x,1.,0.)", vars), + std::runtime_error); + } +} + +namespace { + +// The real-world formula corpus (see corpus.txt in the RooFit JIT-free +// formula planning notes): every RooFormulaVar/RooGenericPdf construction +// string in roofit/*/test/, factory strings, the RooFit tutorials, formulas +// RooFit itself generates (HistFactory, RooProdPdf, convolution bases, ...), +// and the RooFit-realistic subset of test/TFormulaParsingTests.h. +// The category-state entry "cat==cat::c1" appears in its processed form. +const char *const kCorpus[] = { + "floor(x / 2.0) + 1.0", + "floor(x / 2.0)", + "floor(x)", + "floor(x) + floor(y)", + "floor(x) + 1.0", + "TMath::Floor(x) + TMath::Ceil(x) + TMath::Abs(x) + TMath::Tan(x) + TMath::ASin(x / 2.) + TMath::ACos(x / 2.) + " + "TMath::ATan(x) + TMath::Pi() + TMath::E()", + "TMath::PiOver2() + TMath::TanH(x) + TMath::SinH(x) + TMath::Log10(2. + x)", + "2 * x[0] * x[1]", + "@0*0.2e-6 + @1*0.1", + "x + y", + "x + 2.0", + "gauss", + "std::exp(-0.5 * (x*x))", + "x + shift", + "std::pow(x,a)", + "(x-5)*(x-5)*1.2", + "x[0]", + "a1 + x", + "a1 + x + a2 *x*x", + "exp(-2.*x)", + "TMath::Gaus(x, 3, 2)", + "x*x*x+1", + "exp(-0.5*x)", + "TMath::Gaus(x, 5, 0.7)", + "TMath::Gaus(x, 8, 0.8)", + "x[0]==1", // processed form of "cat==cat::c1" + "cat", + "catIndex > 0.5", + "(1+0.1*abs(x)+sin(sqrt(abs(x*alpha+0.1))))", + "sqrt(mean2)", + "a0-a1*sqrt(10*abs(y))", + "0.1*x", + "0.9*x", + "0.0*y", + "0.1*y*y", + "log10(@0)-log10(@1)", + "(x*x+10)", + "x*x+10", + "(1-a)+a*cos((x-c)/b)", + "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))", + "0.5*(std::erf((t-1)/0.5)+1)", + "exp(-@0/ @1)*cosh(@0*@2/2)", + "exp(-@0/ @1)*sinh(@0*@2/2)", + "@1/@0", + "@0*@1*(1-2*@2)", + "@0", + "x - x + 1.0", + "y - y + 1.0", + "0.1 + x*(a + b*x)", + "0.1 + x*(a + x*(b + x*(c + d * x)))", + "log(a*x)", + "ROOT::Math::breitwigner_pdf(x, b, a)", + "ROOT::Math::gaussian_pdf(x, s, m)", + "ROOT::Math::gaussian_pdf(theX, 1, 0)", + "x*std::sqrt(x) + y*std::sqrt(y) + x*y", + "x", + "std::exp(-0.5*(x - mean1)^2/width^2)", + "std::exp(-0.5*(x - mean2)^2/width^2)", + "delta/(sigma*std::sqrt(TMath::TwoPi()))*std::exp(-0.5*(gamma+delta*TMath::ASinH((mass-mu)/sigma))*(gamma+delta*" + "TMath::ASinH((mass-mu)/sigma)))/std::sqrt(1+(mass-mu)*(mass-mu)/(sigma*sigma))", + "var*(par + 1)", + "std::exp(-0.5*(x - mean) * (x - mean) / (sigma * sigma))", + "x * y", + "r + B + y", + "2.7*@0", + "x + 0", + "b*(y<100)", + "1.0 + 1.0*pow(@0,1) + 1.0*pow(@0,2)", + "1.0 + x - x + y - y", + "nbkg_func + 0*x", + "-x[0]", + "x[0] - x[0] + 1", + "1 - x", + "1 + x - x", + "x + y + z", + "mu+shift", + "sigma*1.5", + "sqrt(@0)", + "2 * @0", + "mu*S+B", + "1+1./sqrt(n_off)", + "1+1./sqrt(y)", + "sqrt(n_off)", + "sqrt(y)", + "sqrt(y0)", + "sig+bkg1", + "sig+bkg2", + "1+0.02*alpha_bkg", + "1+0.02*alpha_bkg_A", + "1+0.05*alpha_bkg_B", + "0.5 * pow(1.2, e1)", + "5 * pow(1.3, b1)", + "2*sig*pow(1.2, beta)", + "eff * sig + bkg", + "0.07 * x + 2.0", + "@0*@2+d", + "a * x + c", + "x[1] * x[0] + x[2]", + "@1 * @0 + @2", + "@0 * 2 * @1 + @2", + "f * 3.0/2000. * x * x + (1 - f) / 20.", + "x*x+1", + "1/mean", + "x^4+5*x^3+2*x^2+x+1", + "1+sin(2*@0)", + "acos(cpsi)", + "abs(mean)(in[i + 1]))) { + maxIndex = std::max(maxIndex, in[i + 1] - '0'); + } else if (in[i] == 'x' && i + 1 < n && in[i + 1] == '[' && + (i == 0 || !std::isalnum(static_cast(in[i - 1])))) { + maxIndex = std::max(maxIndex, std::atoi(in.c_str() + i + 2)); + } + } + + std::map named; + std::string out; + std::size_t i = 0; + auto isIdentChar = [](char c) { return std::isalnum(static_cast(c)) || c == '_'; }; + while (i < n) { + const char c = in[i]; + if (std::isdigit(static_cast(c)) || + (c == '.' && i + 1 < n && std::isdigit(static_cast(in[i + 1])))) { + // numeric literal (also hex): copy verbatim + if (c == '0' && i + 1 < n && (in[i + 1] == 'x' || in[i + 1] == 'X')) { + out += in[i++]; + out += in[i++]; + while (i < n && std::isalnum(static_cast(in[i]))) + out += in[i++]; + continue; + } + while (i < n && (std::isdigit(static_cast(in[i])) || in[i] == '.')) + out += in[i++]; + if (i < n && (in[i] == 'e' || in[i] == 'E')) { + std::size_t k = i + 1; + if (k < n && (in[k] == '+' || in[k] == '-')) + ++k; + if (k < n && std::isdigit(static_cast(in[k]))) { + while (i < k) + out += in[i++]; + while (i < n && std::isdigit(static_cast(in[i]))) + out += in[i++]; + } + } + continue; + } + if (c == '@' && i + 1 < n && std::isdigit(static_cast(in[i + 1]))) { + out += "x["; + ++i; + while (i < n && std::isdigit(static_cast(in[i]))) + out += in[i++]; + out += "]"; + continue; + } + if (std::isalpha(static_cast(c)) || c == '_') { + std::size_t start = i; + while (i < n && isIdentChar(in[i])) + ++i; + while (i + 2 < n && in[i] == ':' && in[i + 1] == ':' && + (std::isalpha(static_cast(in[i + 2])) || in[i + 2] == '_')) { + i += 2; + while (i < n && isIdentChar(in[i])) + ++i; + } + const std::string name = in.substr(start, i - start); + if (name == "x" && i < n && in[i] == '[') { + out += name; + while (i < n && in[i] != ']') + out += in[i++]; + if (i < n) + out += in[i++]; // ']' + continue; + } + // function call? + std::size_t j = i; + while (j < n && std::isspace(static_cast(in[j]))) + ++j; + if (j < n && in[j] == '(') { + out += name; + continue; + } + // named variable + auto it = named.find(name); + if (it == named.end()) + it = named.emplace(name, ++maxIndex).first; + out += "x[" + std::to_string(it->second) + "]"; + continue; + } + out += in[i++]; + } + nVars = maxIndex + 1; + return out; +} + +} // namespace + +// First differential smoke test (the full differential campaign is Phase 3): +// every corpus entry that parses on the JIT-free path must evaluate bitwise +// identically to a directly constructed TFormula on random inputs. Also +// report the corpus coverage of the JIT-free path. +TEST(RooFormulaEvaluator, DifferentialCorpus) +{ + int nTotal = 0; + int nAst = 0; + std::vector fallbacks; + + for (const char *entry : kCorpus) { + int nVars = 0; + const std::string processed = normalizeCorpusEntry(entry, nVars); + + TFormula ref("ref", processed.c_str(), /*addToGlobList=*/false); + ASSERT_TRUE(ref.IsValid()) << "corpus entry no longer valid in TFormula: " << entry + << "\n processed: " << processed; + ++nTotal; + + auto prog = RooFormulaParser::compile(processed, nVars); + if (!prog) { + fallbacks.push_back(entry); + continue; + } + ++nAst; + RooExprEvaluator ast{prog}; + + std::mt19937 rng{1234u + static_cast(nTotal)}; + std::uniform_real_distribution dist{-3.0, 3.0}; + for (int trial = 0; trial < 5; ++trial) { + std::vector pars(std::max(nVars, 1)); + for (double &p : pars) { + p = trial == 0 ? 0.5 : trial == 1 ? 2.0 : dist(rng); + } + const double a = ast.eval(pars.data()); + const double t = ref.EvalPar(pars.data()); + EXPECT_TRUE(sameBits(a, t)) << "AST and TFormula disagree for: " << entry << "\n processed: " << processed + << "\n ast = " << std::hexfloat << a << " tformula = " << t << std::defaultfloat; + } + } + + const double coverage = static_cast(nAst) / nTotal; + std::cout << "JIT-free evaluator corpus coverage: " << nAst << "/" << nTotal << " = " << 100. * coverage << "%\n"; + for (auto const &f : fallbacks) { + std::cout << " fallback: " << f << "\n"; + } + EXPECT_GE(coverage, 0.9); +} From d8e80c90ea26e87fd70f09162a87576fe7ddce10 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 27 Aug 2026 16:27:06 +0000 Subject: [PATCH 2/5] [RF] Emit C++ from RooFormula AST for codegen and AD MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Make the JIT-free formula evaluator emit the expression as C++ source, and use that in the codegen backend: the two codegenImpl overloads for RooFormulaVar and RooGenericPdf now inline the emitted expression into the generated code (with the dependents' result names substituted for the formula variables) instead of calling a cling-JIT-compiled TFormula function by name. This removes the per-formula JIT compilation from the codegen path, and Clad now sees plain arithmetic instead of a call across a JIT boundary. The C++ is emitted by walking the same instruction vector that eval() interprets, with explicit parenthesization and with numeric literals formatted at max_digits10 precision, so the emitted code evaluates bitwise identically to the AST evaluator (verified by a new differential test that compiles the emitted code for the corpus expressions and compares against AST evaluation). Formulas that fall back to the TFormula backend keep the previous codegen behavior; the JIT'd function name now flows through a narrow uniqueFuncName() accessor on the evaluator interface, and RooFormulaVar and RooGenericPdf expose which backend handled the formula through a new formulaUsesAstBackend() query (their public getUniqueFuncName() returns an empty string exactly when the expression is inlined instead). getTFormula() is deleted from RooFormula and from the evaluator interface entirely, so no future caller can reintroduce the TFormula coupling, and the temporary lazily-created TFormula stopgap for codegen is gone with it. 🤖 Done with the help of AI --- roofit/codegen/src/CodegenImpl.cxx | 22 ++ roofit/roofitcore/inc/RooFormulaVar.h | 11 +- roofit/roofitcore/inc/RooGenericPdf.h | 11 +- roofit/roofitcore/src/RooExprEvaluator.cxx | 148 +++++++++++++- roofit/roofitcore/src/RooExprEvaluator.h | 10 + roofit/roofitcore/src/RooFormulaEvaluator.h | 29 ++- roofit/roofitcore/src/RooFormulaUtils.cxx | 10 +- roofit/roofitcore/src/RooFormulaVar.cxx | 54 ++--- roofit/roofitcore/src/RooGenericPdf.cxx | 46 +++-- .../roofitcore/src/RooTFormulaEvaluator.cxx | 13 ++ roofit/roofitcore/src/RooTFormulaEvaluator.h | 4 +- .../test/testRooFormulaEvaluator.cxx | 193 +++++++++++++++++- 12 files changed, 485 insertions(+), 66 deletions(-) diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index bc3bdc3556d3e..4c14fcc49a485 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -512,6 +512,20 @@ void codegenImpl(RooGamma &arg, CodegenContext &ctx) void codegenImpl(RooFormulaVar &arg, CodegenContext &ctx) { + // If the formula is handled by RooFit's JIT-free expression backend, the + // expression is inlined into the generated code with the dependents' + // result names substituted for the formula variables. No TFormula (and no + // per-formula JIT compilation) is involved, and Clad sees plain arithmetic + // instead of a call across a JIT boundary. + RooArgList const &deps = arg.dependents(); + std::string expr = arg.emitFormulaCpp([&](unsigned int i) { return ctx.getResult(deps[i]); }); + if (!expr.empty()) { + ctx.addResult(&arg, expr); + return; + } + + // Fallback for formulas the JIT-free backend does not support: call the + // cling-JIT-compiled TFormula function by name. arg.getVal(); // to trigger the creation of the TFormula std::string funcName = arg.getUniqueFuncName(); ctx.collectFunction(funcName); @@ -558,6 +572,14 @@ void codegenImpl(RooGaussian &arg, CodegenContext &ctx) void codegenImpl(RooGenericPdf &arg, CodegenContext &ctx) { + // See the comments in codegenImpl(RooFormulaVar&, CodegenContext&). + RooArgList const &deps = arg.dependents(); + std::string expr = arg.emitFormulaCpp([&](unsigned int i) { return ctx.getResult(deps[i]); }); + if (!expr.empty()) { + ctx.addResult(&arg, expr); + return; + } + arg.getVal(); // to trigger the creation of the TFormula std::string funcName = arg.getUniqueFuncName(); ctx.collectFunction(funcName); diff --git a/roofit/roofitcore/inc/RooFormulaVar.h b/roofit/roofitcore/inc/RooFormulaVar.h index b7cf0232090a3..a6a811b82d8f6 100644 --- a/roofit/roofitcore/inc/RooFormulaVar.h +++ b/roofit/roofitcore/inc/RooFormulaVar.h @@ -21,6 +21,7 @@ #include "RooListProxy.h" #include "RooAbsBinning.h" +#include #include #include #include @@ -29,7 +30,6 @@ class RooArgSet ; class RooFormulaEvaluator; class RooAbsRealLValue; -class TFormula; class RooFormulaVar : public RooAbsReal { public: @@ -84,7 +84,15 @@ class RooFormulaVar : public RooAbsReal { double evaluate() const override ; void doEval(RooFit::EvalContext &ctx) const override; + /// Name of the cling-JIT-compiled function that evaluates this formula, + /// which generated code from the codegen fallback path calls by name. + /// \note Returns an empty string when the formula is handled by the JIT-free + /// formula backend (the default for supported expressions, see + /// formulaUsesAstBackend()); codegen then inlines the expression via + /// emitFormulaCpp() instead. std::string getUniqueFuncName() const; + std::string emitFormulaCpp(std::function const &varName) const; + bool formulaUsesAstBackend() const; std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const override; @@ -97,7 +105,6 @@ class RooFormulaVar : public RooAbsReal { RooListProxy _actualVars; ///< Actual parameters used by formula engine mutable std::unique_ptr _evaluator; /// _tFormulaForCodegen; /// #include #include #include class RooArgList ; class RooFormulaEvaluator; -class TFormula; class RooAbsRealLValue; class RooGenericPdf : public RooAbsPdf { @@ -66,7 +66,15 @@ class RooGenericPdf : public RooAbsPdf { const char* expression() const { return _formExpr.Data(); } const RooArgList& dependents() const { return _actualVars; } + /// Name of the cling-JIT-compiled function that evaluates this formula, + /// which generated code from the codegen fallback path calls by name. + /// \note Returns an empty string when the formula is handled by the JIT-free + /// formula backend (the default for supported expressions, see + /// formulaUsesAstBackend()); codegen then inlines the expression via + /// emitFormulaCpp() instead. std::string getUniqueFuncName() const; + std::string emitFormulaCpp(std::function const &varName) const; + bool formulaUsesAstBackend() const; void setBinning(const RooAbsRealLValue &obs, const RooAbsBinning &binning, bool checkFlatness = true); const RooAbsBinning *getBinning(const RooAbsRealLValue &obs) const; @@ -87,7 +95,6 @@ class RooGenericPdf : public RooAbsPdf { bool isValidReal(double /*value*/, bool /*printError*/) const override { return true; } mutable std::unique_ptr _evaluator; /// _tFormulaForCodegen; ///> _binnings; ///< User-defined binnings, keyed by the observable's index diff --git a/roofit/roofitcore/src/RooExprEvaluator.cxx b/roofit/roofitcore/src/RooExprEvaluator.cxx index 51695d0bc5554..61d075c4b397f 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.cxx +++ b/roofit/roofitcore/src/RooExprEvaluator.cxx @@ -16,6 +16,10 @@ #include #include +#include +#include +#include +#include #include #include @@ -33,20 +37,23 @@ Entry F0(const char *name, double (*fn)()) return e; } -Entry F1(const char *name, double (*fn)(double), TypeRule rule = TypeRule::Double) +Entry F1(const char *name, double (*fn)(double), TypeRule rule = TypeRule::Double, const char *cppName = nullptr) { Entry e; e.name = name; + e.cppName = cppName; e.arity = 1; e.rule = rule; e.fn1 = fn; return e; } -Entry F2(const char *name, double (*fn)(double, double), TypeRule rule = TypeRule::Double) +Entry F2(const char *name, double (*fn)(double, double), TypeRule rule = TypeRule::Double, + const char *cppName = nullptr) { Entry e; e.name = name; + e.cppName = cppName; e.arity = 2; e.rule = rule; e.fn2 = fn; @@ -154,8 +161,10 @@ std::vector makeTable() F1("std::abs", +[](double x) { return std::fabs(x); }, TypeRule::SameAsFirstArg), F1("fabs", +[](double x) { return std::fabs(x); }), F1("std::fabs", +[](double x) { return std::fabs(x); }), - F1("int", castInt, TypeRule::Int), // C++ functional cast: truncation towards zero - F1("sq", square), // TFormula shortcut for TMath::Sq(Double_t) + // C++ functional cast: truncation towards zero + F1("int", castInt, TypeRule::Int, "int"), + // TFormula shortcut for TMath::Sq(Double_t) + F1("sq", square, TypeRule::Double, "TMath::Sq"), // one-argument functions, TMath spellings F1("TMath::Sqrt", +[](double x) { return TMath::Sqrt(x); }), F1("TMath::Exp", +[](double x) { return TMath::Exp(x); }), @@ -195,7 +204,8 @@ std::vector makeTable() F2("max", stdMax, TypeRule::MinMax), F2("std::max", stdMax, TypeRule::MinMax), F2("TMath::Max", tmathMax, TypeRule::MinMax), - F2("sign", sign, TypeRule::SameAsFirstArg), // TFormula shortcut for TMath::Sign + // TFormula shortcut for TMath::Sign + F2("sign", sign, TypeRule::SameAsFirstArg, "TMath::Sign"), F2("TMath::Sign", sign, TypeRule::SameAsFirstArg), // zero-argument constants (folded to Op::Const at parse time) F0("TMath::Pi", +[]() { return TMath::Pi(); }), @@ -342,4 +352,132 @@ double RooExprEvaluator::eval(const double *vars) const return stack[0]; } +namespace { + +/// Format a double as a C++ expression of type double that parses back to the +/// exact same value: max_digits10 (17) significant decimal digits guarantee +/// the bitwise round-trip (this is also the formatting convention used +/// elsewhere in RooFit codegen). +std::string emitDouble(double val) +{ + if (std::isnan(val)) { + return "std::numeric_limits::quiet_NaN()"; + } + if (std::isinf(val)) { + return val > 0 ? "std::numeric_limits::infinity()" : "(-std::numeric_limits::infinity())"; + } + std::stringstream ss; + // The formatting must not depend on the global locale: a comma decimal + // separator (e.g. from a German locale) would corrupt the emitted C++. + ss.imbue(std::locale::classic()); + ss << std::setprecision(std::numeric_limits::max_digits10) << val; + std::string out = ss.str(); + // The emitted literal must have type double: an integer-looking literal + // like `2` would be an int in C++, with different division semantics. + if (out.find_first_of(".eE") == std::string::npos) { + out += ".0"; + } + // Parse-time constant folding (e.g. of TMath::Pi()) can in principle + // produce negative values; keep every stack entry self-contained. + if (out[0] == '-') { + out = "(" + out + ")"; + } + return out; +} + +/// The C++ spelling emitted for a function-table entry (see Entry::cppName). +std::string emissionName(RooFormulaFunctions::Entry const &entry) +{ + if (entry.cppName) { + return entry.cppName; + } + const std::string name = entry.name; + return name.find("::") == std::string::npos ? "std::" + name : name; +} + +} // namespace + +//////////////////////////////////////////////////////////////////////////////// +/// Emit the expression as C++ source by walking the same instruction sequence +/// that eval() interprets, so evaluation and code generation cannot diverge +/// structurally. Every intermediate result is kept fully parenthesized, so no +/// operator-precedence reasoning is involved. The emitted operators and +/// function spellings reproduce what the cling-JIT-compiled TFormula code +/// called for the same formula, so all three semantic paths (AST evaluation, +/// emitted C++, TFormula fallback) agree bitwise. +std::string RooExprEvaluator::emitCpp(std::function const &varName) const +{ + std::vector stack; + auto const *funcs = RooFormulaFunctions::table(); + + auto pop = [&]() { + std::string out = std::move(stack.back()); + stack.pop_back(); + return out; + }; + auto binary = [&](const char *sym) { + const std::string b = pop(); + stack.back() = "(" + stack.back() + " " + sym + " " + b + ")"; + }; + + for (Instr const &ins : _program->code) { + switch (ins.op) { + case Op::Const: stack.push_back(emitDouble(ins.konst)); break; + case Op::Var: stack.push_back("(" + varName(ins.arg) + ")"); break; + case Op::Add: binary("+"); break; + case Op::Sub: binary("-"); break; + case Op::Mul: binary("*"); break; + case Op::Div: binary("/"); break; + case Op::Neg: stack.back() = "(-" + stack.back() + ")"; break; + case Op::Not: stack.back() = "(!" + stack.back() + ")"; break; + case Op::LT: binary("<"); break; + case Op::LE: binary("<="); break; + case Op::GT: binary(">"); break; + case Op::GE: binary(">="); break; + case Op::EQ: binary("=="); break; + case Op::NE: binary("!="); break; + // In C++, `&&`/`||` short-circuit and `?:` evaluates only the taken + // branch, while eval() always evaluates both operands. The resulting + // values are identical; the difference is only observable if + // floating-point exceptions are trapped. + case Op::And: binary("&&"); break; + case Op::Or: binary("||"); break; + case Op::Select: { + const std::string b = pop(); + const std::string a = pop(); + stack.back() = "(" + stack.back() + " ? " + a + " : " + b + ")"; + break; + } + case Op::Pow: { + const std::string b = pop(); + stack.back() = "std::pow(" + stack.back() + ", " + b + ")"; + break; + } + case Op::Sq: stack.back() = "TMath::Sq(" + stack.back() + ")"; break; + case Op::IntNorm: stack.back() = "(" + stack.back() + " + 0.0)"; break; + case Op::Call1: stack.back() = emissionName(funcs[ins.arg]) + "(" + stack.back() + ")"; break; + case Op::Call2: { + const std::string b = pop(); + stack.back() = emissionName(funcs[ins.arg]) + "(" + stack.back() + ", " + b + ")"; + break; + } + case Op::Call3: { + const std::string c = pop(); + const std::string b = pop(); + stack.back() = emissionName(funcs[ins.arg]) + "(" + stack.back() + ", " + b + ", " + c + ")"; + break; + } + case Op::Call4: { + const std::string d = pop(); + const std::string c = pop(); + const std::string b = pop(); + stack.back() = emissionName(funcs[ins.arg]) + "(" + stack.back() + ", " + b + ", " + c + ", " + d + ")"; + break; + } + } + } + + return stack.back(); +} + /// \endcond diff --git a/roofit/roofitcore/src/RooExprEvaluator.h b/roofit/roofitcore/src/RooExprEvaluator.h index 45d743199940b..b16ee9ec3c29c 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.h +++ b/roofit/roofitcore/src/RooExprEvaluator.h @@ -39,6 +39,12 @@ enum class TypeRule : std::uint8_t { struct Entry { const char *name = nullptr; ///< accepted spelling in the formula + /// Spelling emitted in generated C++ for this entry. A nullptr means + /// "derive from name": qualified names (TMath::Erf, std::sin) are emitted + /// as-is, bare libm/std names get a std:: qualification (sin -> std::sin, + /// resolving to the same function the JIT-compiled code called). Only + /// entries whose emission cannot be derived this way set it explicitly. + const char *cppName = nullptr; std::uint8_t arity = 0; TypeRule rule = TypeRule::Double; double (*fn0)() = nullptr; @@ -123,6 +129,10 @@ class RooExprEvaluator final : public RooFormulaEvaluator { /// The processed formula string this program was compiled from. std::string processedFormula() const { return _program->formula; } + bool canEmitCpp() const override { return true; } + + std::string emitCpp(std::function const &varName) const override; + private: std::shared_ptr _program; }; diff --git a/roofit/roofitcore/src/RooFormulaEvaluator.h b/roofit/roofitcore/src/RooFormulaEvaluator.h index fe922e62c22fe..4fbe578773ffb 100644 --- a/roofit/roofitcore/src/RooFormulaEvaluator.h +++ b/roofit/roofitcore/src/RooFormulaEvaluator.h @@ -13,9 +13,9 @@ #ifndef ROO_FORMULA_EVALUATOR #define ROO_FORMULA_EVALUATOR +#include #include - -class TFormula; +#include /// Abstract interface for evaluating a processed formula expression, i.e. one /// normalized by RooFormulaUtils::processFormula() to the `x[i]`-only dialect, @@ -30,11 +30,26 @@ class RooFormulaEvaluator { /// Return a deep copy of this evaluator. virtual std::unique_ptr clone() const = 0; - /// Return the underlying TFormula. Only the TFormula-backed evaluator - /// returns a non-nullptr. This accessor only exists to support the - /// getUniqueFuncName() functions used by the codegen backend and will be - /// removed together with them. - virtual TFormula *getTFormula() const { return nullptr; } + /// Whether emitCpp() can emit this expression as C++ source. Only the + /// JIT-free expression backend can; the TFormula backend cannot. + virtual bool canEmitCpp() const { return false; } + + /// Emit this expression as C++ source with explicit parenthesization, for + /// RooFit code generation and automatic differentiation. `varName(i)` + /// supplies the emitted name for `x[i]`. Returns an empty string if this + /// evaluator cannot emit C++ (see canEmitCpp()); the codegen caller then + /// uses the TFormula fallback path via uniqueFuncName(). + virtual std::string emitCpp(std::function const & /*varName*/) const { return {}; } + + /// Name of the cling-JIT-compiled function that evaluates this formula. + /// Only meaningfully implemented by the TFormula backend, where it serves + /// the codegen fallback path for formulas that cannot emitCpp(): the + /// generated code calls that function by name. Empty otherwise. + virtual std::string uniqueFuncName() const { return {}; } + + /// Propagate a rename of the owning object to any named objects held + /// by the evaluator (the TFormula backend renames its TFormula). + virtual void setName(const char * /*name*/) {} }; #endif diff --git a/roofit/roofitcore/src/RooFormulaUtils.cxx b/roofit/roofitcore/src/RooFormulaUtils.cxx index 461910968aea7..0e156e851f2ba 100644 --- a/roofit/roofitcore/src/RooFormulaUtils.cxx +++ b/roofit/roofitcore/src/RooFormulaUtils.cxx @@ -48,8 +48,6 @@ unsupported expressions into hard errors (useful for testing). #include "RooFormulaParser.h" #include "RooTFormulaEvaluator.h" -#include "TFormula.h" - #include #include #include @@ -501,15 +499,13 @@ RooFormulaEvaluator &RooFormulaUtils::ensureEvaluator(std::unique_ptr RooFormulaUtils::cloneEvaluator(RooFormulaEvaluator const &other, const char *newName) { std::unique_ptr out = other.clone(); - if (TFormula *tFormula = out->getTFormula()) { - tFormula->SetName(newName); - } + out->setName(newName); return out; } diff --git a/roofit/roofitcore/src/RooFormulaVar.cxx b/roofit/roofitcore/src/RooFormulaVar.cxx index c19c1c2e019c4..9e488e92bcd01 100644 --- a/roofit/roofitcore/src/RooFormulaVar.cxx +++ b/roofit/roofitcore/src/RooFormulaVar.cxx @@ -51,8 +51,6 @@ #include "RooFormulaUtils.h" #include "RooAbsRealLValue.h" -#include "TFormula.h" - #ifdef ROOFIT_LEGACY_EVAL_BACKEND #include "RooNLLVar.h" #include "RooChi2Var.h" @@ -327,29 +325,39 @@ double RooFormulaVar::defaultErrorLevel() const } //////////////////////////////////////////////////////////////////////////////// -/// Name of the cling-JIT-compiled function that evaluates this formula, which -/// generated code from the codegen path calls by name. -/// -/// If the formula is evaluated by the TFormula backend, this is the function -/// of the evaluating TFormula itself. Otherwise, a TFormula is created lazily, -/// only to serve the codegen path. -/// -/// TODO(Phase 2.5): remove the lazily-created TFormula once codegen emits C++ -/// for the expression directly from the parsed representation instead of -/// calling the JIT-compiled TFormula function. Until then, codegen keeps -/// working exactly as before, at the cost of one JIT compilation per formula -/// -- but only when codegen is actually used. Like codegen itself, this lazy -/// creation is not thread-safe. +/// Name of the cling-JIT-compiled function evaluating this formula, which the +/// codegen fallback path calls by name in generated code. Empty when the +/// formula is evaluated by the JIT-free expression backend (codegen then +/// inlines the expression via emitFormulaCpp() instead). std::string RooFormulaVar::getUniqueFuncName() const { - if (TFormula *tFormula = evaluator().getTFormula()) { - return tFormula->GetUniqueFuncName().Data(); - } - if (!_tFormulaForCodegen) { - // evaluator() above has normalized _formExpr to the processed `x[i]` dialect. - _tFormulaForCodegen = std::make_unique(GetName(), _formExpr.Data(), /*addToGlobList=*/false); - } - return _tFormulaForCodegen->GetUniqueFuncName().Data(); + return evaluator().uniqueFuncName(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// If the formula expression can be emitted as inline C++ (i.e. it is +/// evaluated by the JIT-free expression backend), return the emitted +/// expression, with `varName(i)` supplying the generated name of +/// `dependents()[i]`. Return an empty string otherwise; codegen then falls +/// back to calling the cling-JIT-compiled TFormula function by name (see +/// getUniqueFuncName()). +std::string RooFormulaVar::emitFormulaCpp(std::function const &varName) const +{ + return evaluator().emitCpp(varName); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Whether the formula expression is evaluated by RooFit's built-in JIT-free +/// (AST) formula backend, which is the default for supported expressions. +/// Returns false if it is evaluated by the TFormula (cling JIT) fallback +/// backend instead, either because the built-in parser does not support the +/// expression or because the TFormula backend was forced with the +/// ROOFIT_FORMULA_BACKEND environment variable. Exactly the formulas on the +/// JIT-free backend can be emitted as inline C++ (emitFormulaCpp()); the +/// others report the name of their JIT-compiled function (getUniqueFuncName()). +bool RooFormulaVar::formulaUsesAstBackend() const +{ + return evaluator().canEmitCpp(); } std::unique_ptr diff --git a/roofit/roofitcore/src/RooGenericPdf.cxx b/roofit/roofitcore/src/RooGenericPdf.cxx index ea411c6706a3f..dd938f20a2788 100644 --- a/roofit/roofitcore/src/RooGenericPdf.cxx +++ b/roofit/roofitcore/src/RooGenericPdf.cxx @@ -35,8 +35,6 @@ class documentation. #include "RooFormulaUtils.h" #include "RooAbsRealLValue.h" -#include "TFormula.h" - using std::istream, std::ostream, std::endl; @@ -228,19 +226,37 @@ void RooGenericPdf::writeToStream(ostream& os, bool compact) const } //////////////////////////////////////////////////////////////////////////////// -/// Name of the cling-JIT-compiled function that evaluates this formula, which -/// generated code from the codegen path calls by name. -/// -/// See RooFormulaVar::getUniqueFuncName() for the details of the lazily -/// created TFormula on the JIT-free expression backend. +/// Name of the cling-JIT-compiled function evaluating this formula, which the +/// codegen fallback path calls by name in generated code. Empty when the +/// formula is evaluated by the JIT-free expression backend (codegen then +/// inlines the expression via emitFormulaCpp() instead). std::string RooGenericPdf::getUniqueFuncName() const { - if (TFormula *tFormula = evaluator().getTFormula()) { - return tFormula->GetUniqueFuncName().Data(); - } - if (!_tFormulaForCodegen) { - // evaluator() above has normalized _formExpr to the processed `x[i]` dialect. - _tFormulaForCodegen = std::make_unique(GetName(), _formExpr.Data(), /*addToGlobList=*/false); - } - return _tFormulaForCodegen->GetUniqueFuncName().Data(); + return evaluator().uniqueFuncName(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// If the formula expression can be emitted as inline C++ (i.e. it is +/// evaluated by the JIT-free expression backend), return the emitted +/// expression, with `varName(i)` supplying the generated name of +/// `dependents()[i]`. Return an empty string otherwise; codegen then falls +/// back to calling the cling-JIT-compiled TFormula function by name (see +/// getUniqueFuncName()). +std::string RooGenericPdf::emitFormulaCpp(std::function const &varName) const +{ + return evaluator().emitCpp(varName); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Whether the formula expression is evaluated by RooFit's built-in JIT-free +/// (AST) formula backend, which is the default for supported expressions. +/// Returns false if it is evaluated by the TFormula (cling JIT) fallback +/// backend instead, either because the built-in parser does not support the +/// expression or because the TFormula backend was forced with the +/// ROOFIT_FORMULA_BACKEND environment variable. Exactly the formulas on the +/// JIT-free backend can be emitted as inline C++ (emitFormulaCpp()); the +/// others report the name of their JIT-compiled function (getUniqueFuncName()). +bool RooGenericPdf::formulaUsesAstBackend() const +{ + return evaluator().canEmitCpp(); } diff --git a/roofit/roofitcore/src/RooTFormulaEvaluator.cxx b/roofit/roofitcore/src/RooTFormulaEvaluator.cxx index d9a72a3df2160..a2068ec4d0bcc 100644 --- a/roofit/roofitcore/src/RooTFormulaEvaluator.cxx +++ b/roofit/roofitcore/src/RooTFormulaEvaluator.cxx @@ -88,4 +88,17 @@ std::unique_ptr RooTFormulaEvaluator::clone() const return std::make_unique(*this); } +//////////////////////////////////////////////////////////////////////////////// +/// Return the name of the cling-JIT-compiled function that evaluates this +/// formula, which the codegen fallback path calls by name in generated code. +std::string RooTFormulaEvaluator::uniqueFuncName() const +{ + return _tFormula->GetUniqueFuncName().Data(); +} + +void RooTFormulaEvaluator::setName(const char *name) +{ + _tFormula->SetName(name); +} + /// \endcond diff --git a/roofit/roofitcore/src/RooTFormulaEvaluator.h b/roofit/roofitcore/src/RooTFormulaEvaluator.h index bb4a1301c1f60..527fce92d57d4 100644 --- a/roofit/roofitcore/src/RooTFormulaEvaluator.h +++ b/roofit/roofitcore/src/RooTFormulaEvaluator.h @@ -35,8 +35,8 @@ class RooTFormulaEvaluator : public RooFormulaEvaluator { double eval(const double *vars) const override; std::unique_ptr clone() const override; - - TFormula *getTFormula() const override { return _tFormula.get(); } + std::string uniqueFuncName() const override; + void setName(const char *name) override; private: std::unique_ptr _tFormula; ///< The formula used to compute values diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx index 0b4cbc27e595b..682d0eb8d9c5c 100644 --- a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -8,11 +8,13 @@ #include "../src/RooExprEvaluator.h" #include +#include #include #include #include #include +#include #include #include @@ -20,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -492,9 +495,31 @@ TEST(RooFormulaEvaluator, RooFormulaIntegration) EXPECT_STREQ(fVar.expression(), "x[0]*x[1]+sin(x[0])"); EXPECT_TRUE(sameBits(fVar.getVal(), 2.0 * 3.0 + std::sin(2.0))); - // the codegen accessor lazily provides a JIT-compiled TFormula function - // even on the AST path - EXPECT_FALSE(fVar.getUniqueFuncName().empty()); + // on the AST path, the expression is emitted as C++ for codegen instead of + // JIT-compiling a TFormula; there is no JIT'd function name + ASSERT_TRUE(f->canEmitCpp()); + auto varName = [](unsigned int i) { return "v[" + std::to_string(i) + "]"; }; + EXPECT_EQ(f->emitCpp(varName), "(((v[0]) * (v[1])) + std::sin((v[0])))"); + EXPECT_TRUE(f->uniqueFuncName().empty()); + EXPECT_TRUE(fVar.getUniqueFuncName().empty()); +} + +// Codegen indexes the pruned list of actually-used dependents. Since the +// owning classes prune unused variables and reindex the expression before the +// evaluation engine is created, emitFormulaCpp() must name the variables +// consistently with dependents(). +TEST(RooFormulaEvaluator, EmitCppDependentRemap) +{ + ScopedBackendEnv env{nullptr}; + + RooRealVar x("x", "x", 2.0); + RooRealVar y("y", "y", 3.0); + + // x (list index 0) is unused, so y is dependents()[0]. + RooFormulaVar f("f", "y*2", {x, y}); + ASSERT_EQ(f.dependents().size(), 1u); + std::string expr = f.emitFormulaCpp([](unsigned int i) { return "dep" + std::to_string(i); }); + EXPECT_EQ(expr, "((dep0) * 2.0)"); } TEST(RooFormulaEvaluator, BackendOverride) @@ -511,6 +536,11 @@ TEST(RooFormulaEvaluator, BackendOverride) auto g = RooFormulaUtils::makeFormulaEvaluator("g", "ROOT::Math::normal_pdf(x,1.,0.)", vars); EXPECT_FALSE(isAstBackend(*g)); EXPECT_TRUE(sameBits(RooFormulaUtils::evalFormula(*g, vars), ROOT::Math::normal_pdf(2.0, 1.0, 0.0))); + // the fallback backend cannot emit C++; codegen instead calls the + // cling-JIT-compiled TFormula function by its unique name + EXPECT_FALSE(g->canEmitCpp()); + EXPECT_TRUE(g->emitCpp([](unsigned int) { return std::string{"v"}; }).empty()); + EXPECT_FALSE(g->uniqueFuncName().empty()); } { ScopedBackendEnv env{"tformula"}; @@ -529,6 +559,35 @@ TEST(RooFormulaEvaluator, BackendOverride) } } +// The public backend query on RooFormulaVar and RooGenericPdf: true when the +// formula is evaluated by the JIT-free AST backend (the default for supported +// expressions), false on the TFormula fallback backend. +TEST(RooFormulaEvaluator, PublicBackendQuery) +{ + RooRealVar a("a", "a", 1.1); + RooRealVar x("x", "x", 2.0); + RooRealVar b("b", "b", 0.3); + + { + ScopedBackendEnv env{nullptr}; + RooFormulaVar f("f", "a*x+b", {a, x, b}); + EXPECT_TRUE(f.formulaUsesAstBackend()); + EXPECT_TRUE(f.getUniqueFuncName().empty()); + RooGenericPdf p("p", "a*x+b", {a, x, b}); + EXPECT_TRUE(p.formulaUsesAstBackend()); + EXPECT_TRUE(p.getUniqueFuncName().empty()); + } + { + ScopedBackendEnv env{"tformula"}; + RooFormulaVar f("f", "a*x+b", {a, x, b}); + EXPECT_FALSE(f.formulaUsesAstBackend()); + EXPECT_FALSE(f.getUniqueFuncName().empty()); + RooGenericPdf p("p", "a*x+b", {a, x, b}); + EXPECT_FALSE(p.formulaUsesAstBackend()); + EXPECT_FALSE(p.getUniqueFuncName().empty()); + } +} + namespace { // The real-world formula corpus (see corpus.txt in the RooFit JIT-free @@ -846,3 +905,131 @@ TEST(RooFormulaEvaluator, DifferentialCorpus) } EXPECT_GE(coverage, 0.9); } + +namespace { + +using EmittedFunc = double (*)(double const *); + +/// Declare the emitted expression as a function of the variable array `v` in +/// the interpreter and return a pointer to the compiled function. +EmittedFunc compileEmitted(std::string const &expr) +{ + static bool headersDeclared = + gInterpreter->Declare("#include \"TMath.h\"\n#include \n#include \n#include \n"); + if (!headersDeclared) { + return nullptr; + } + static int counter = 0; + const std::string fname = "rooFormulaEmitTestFunc" + std::to_string(counter++); + const std::string code = "double " + fname + "(double const *v) { return " + expr + "; }"; + if (!gInterpreter->Declare(code.c_str())) { + return nullptr; + } + return reinterpret_cast(gInterpreter->Calc(("(void *) " + fname).c_str())); +} + +} // namespace + +// Emitted-code agreement (Phase 3 item 5, brought forward minimally): for +// every corpus expression the JIT-free parser accepts, emit the C++ +// expression, compile it with the interpreter, and require bitwise agreement +// with AST evaluation on random inputs. The C++ is emitted from the same +// instruction vector that eval() walks, so this directly validates the +// emission itself: the operator spellings, the function-name mapping table, +// and the exact round-trip of numeric literals. +TEST(RooFormulaEvaluator, EmittedCppDifferential) +{ + auto varName = [](unsigned int i) { return "v[" + std::to_string(i) + "]"; }; + + int iEntry = 0; + for (const char *entry : kCorpus) { + ++iEntry; + int nVars = 0; + const std::string processed = normalizeCorpusEntry(entry, nVars); + + auto prog = RooFormulaParser::compile(processed, nVars); + if (!prog) { + continue; // fallback expressions have no emission; covered elsewhere + } + RooExprEvaluator ast{prog}; + + const std::string expr = ast.emitCpp(varName); + ASSERT_FALSE(expr.empty()) << entry; + EmittedFunc fn = compileEmitted(expr); + ASSERT_NE(fn, nullptr) << "emitted C++ failed to compile for: " << entry << "\n emitted: " << expr; + + std::mt19937 rng{987u + static_cast(iEntry)}; + std::uniform_real_distribution dist{-3.0, 3.0}; + for (int trial = 0; trial < 5; ++trial) { + std::vector pars(std::max(nVars, 1)); + for (double &p : pars) { + p = trial == 0 ? 0.5 : trial == 1 ? 2.0 : dist(rng); + } + const double a = ast.eval(pars.data()); + const double c = fn(pars.data()); + EXPECT_TRUE(sameBits(a, c)) << "AST and emitted C++ disagree for: " << entry << "\n emitted: " << expr + << "\n ast = " << std::hexfloat << a << " emitted = " << c << std::defaultfloat; + } + } +} + +// Numeric literals must survive emit -> compile -> eval bitwise, so they are +// emitted with max_digits10 (17) significant digits. A lossy emission (e.g. +// the default 6-digit %g formatting) would show up here: 0.1 and friends are +// not exactly representable. +TEST(RooFormulaEvaluator, EmittedLiteralRoundTrip) +{ + auto varName = [](unsigned int) { return std::string{"v[0]"}; }; + + auto prog = RooFormulaParser::compile("x[0]*0.1+0.2e-6*3.360779", 1); + ASSERT_TRUE(prog); + RooExprEvaluator ast{prog}; + const std::string expr = ast.emitCpp(varName); + // 0.1 must be emitted with enough digits for an exact round-trip + EXPECT_NE(expr.find("0.10000000000000001"), std::string::npos) << expr; + EmittedFunc fn = compileEmitted(expr); + ASSERT_NE(fn, nullptr) << expr; + for (double v : {0.3, 1.0, 7.7, 1e30, 1e-30, -2.5}) { + EXPECT_TRUE(sameBits(ast.eval(&v), fn(&v))) << expr << " at v = " << v; + } + + // integer-spelled literals must be emitted with double type: `2` in the + // formula dialect is emitted as `2.0` + auto prog2 = RooFormulaParser::compile("7./2", 0); + ASSERT_TRUE(prog2); + EXPECT_EQ(RooExprEvaluator{prog2}.emitCpp(varName), "(7.0 / 2.0)"); +} + +// The emitted C++ must not depend on the global locale: under a comma-decimal +// locale a default-constructed stream would format 0.5 as "0,5", corrupting +// the generated code. No comma-decimal OS locale is installed on every test +// machine, so the locale is built from a custom numpunct facet instead. The +// global locale is restored before any assertion can bail out of the test. +TEST(RooFormulaEvaluator, EmitCppLocaleIndependent) +{ + struct CommaPunct : std::numpunct { + char do_decimal_point() const override { return ','; } + }; + + const std::locale old = std::locale::global(std::locale{std::locale::classic(), new CommaPunct}); + std::string expr; + std::string formatted; + try { + std::stringstream ss; // sanity check: the facet is actually in effect + ss << 0.5; + formatted = ss.str(); + auto prog = RooFormulaParser::compile("x[0]*0.5+1.25", 1); + if (prog) { + expr = RooExprEvaluator{prog}.emitCpp([](unsigned int i) { return "v[" + std::to_string(i) + "]"; }); + } + } catch (...) { + std::locale::global(old); + throw; + } + std::locale::global(old); + + EXPECT_EQ(formatted, "0,5"); + EXPECT_NE(expr.find("0.5"), std::string::npos) << expr; + EXPECT_NE(expr.find("1.25"), std::string::npos) << expr; + EXPECT_EQ(expr.find(','), std::string::npos) << expr; +} From 0e79371fb51e1001c845f45896a782158d388b8c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 28 Aug 2026 07:26:51 +0000 Subject: [PATCH 3/5] [RF] Add differential and workspace-IO tests for the formula AST evaluator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 3 (differential testing) and phase 3.5 (workspace I/O) of the JIT-free RooFormula evaluation work: * A deterministic random-expression generator over the full supported grammar (all operators including `^`/`**`, comparisons, logical operators, ternary, every allow-listed function spelling sampled from the actual table, and edge-case inputs including 0, +-1, 1e300, 1e-300 and negative arguments to sqrt/log). 500 expressions x 8 input vectors must agree bitwise with TFormula/cling, and a 150-expression sample additionally goes through emitCpp() + interpreter compilation and must agree bitwise with AST evaluation. No tolerances anywhere. * The generator found four real divergences between the AST dialect and what TFormula/cling actually compiles; the parser now falls back to the TFormula backend for all of them, keeping today's behavior: - a textually adjacent `++` is TFormula's linear-combination separator, not an addition (and runs of three or more `-` survive TFormula's double-negation rewrite as a pre-decrement that does not compile); - bare chained comparisons (`a < b < c`) are invalid in TFormula because cling compiles with -Wparentheses promoted to an error; - `^` with an explicit sign on a parenthesized exponent is broken in TFormula (`x^-(a+b)` compiles as pow(x,-(a)+b)); - sign()/TMath::Sign with a bool-typed first argument resolves to the generic TMath::Sign template returning bool, which is not copysign. Expression typing is now tracked as double/int/bool to detect this (and bool/int min/max mixes, which do not compile). * The same typing information also folds int-typed constant subexpressions in int64 at parse time, because cling computes them in wrapping int32 arithmetic ("100000*100000" is 1410065408 there, not 1e10). Any int-typed constant intermediate that leaves the int32 range falls back, as do integer literals too large for int (which have type long or unsigned int in C++, not int). * testRooFuncWrapper gains a gradient-agreement test: for a Gaussian with RooFormulaVar parameters and for a numerically-integrated RooGenericPdf, the Clad gradient from the new inlined-expression codegen path must match the gradient from the old call-the-JIT'd-TFormula path to 1e-9 and numerical differentiation to 1e-4. * New testRooFormulaEvaluatorIO.cxx covers workspace persistence: - round-trip through the AST path (the direct regression test for the formulaString()/_formExpr persistence landmine), asserting non-empty persisted strings and bitwise-identical evaluation after reading back; - reading a legacy workspace fixture written by a pre-AST build (testRooFormulaEvaluator_legacy_ws.root), asserting bitwise agreement with the recorded 17-digit reference values on both backends; - on-disk fidelity: writing the same content through the AST and the TFormula backends produces identical persisted formula strings and identical streamed class versions, both matching the legacy file. The real-world corpus coverage of the AST path is unchanged at 152/156 (97.4%), reported via GTest properties with a hard floor of 90%. 🤖 Done with the help of AI --- roofit/roofitcore/src/RooExprEvaluator.cxx | 13 +- roofit/roofitcore/src/RooExprEvaluator.h | 16 +- roofit/roofitcore/src/RooFormulaParser.cxx | 218 ++++++-- roofit/roofitcore/test/CMakeLists.txt | 4 +- .../test/testRooFormulaEvaluator.cxx | 496 ++++++++++++++++++ .../test/testRooFormulaEvaluatorIO.cxx | 296 +++++++++++ .../testRooFormulaEvaluator_legacy_ws.root | Bin 0 -> 9243 bytes roofit/roofitcore/test/testRooFuncWrapper.cxx | 116 ++++ 8 files changed, 1114 insertions(+), 45 deletions(-) create mode 100644 roofit/roofitcore/test/testRooFormulaEvaluatorIO.cxx create mode 100644 roofit/roofitcore/test/testRooFormulaEvaluator_legacy_ws.root diff --git a/roofit/roofitcore/src/RooExprEvaluator.cxx b/roofit/roofitcore/src/RooExprEvaluator.cxx index 61d075c4b397f..8044aca89f345 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.cxx +++ b/roofit/roofitcore/src/RooExprEvaluator.cxx @@ -104,8 +104,11 @@ std::vector makeTable() auto tmathMin = +[](double a, double b) { return TMath::Min(a, b); }; auto tmathMax = +[](double a, double b) { return TMath::Max(a, b); }; // Both the `sign` shortcut and TMath::Sign resolve to TMath::Sign, which is - // std::copysign for double arguments. (For an integer first argument cling - // picks the TMath::Sign template, whose result is numerically identical.) + // std::copysign for double arguments. (For an int first argument cling picks + // the generic TMath::Sign template, which also follows the sign bit of the + // second argument and is numerically identical for int values. A *bool* + // first argument is rejected by the parser: the template would return bool, + // turning e.g. Sign(true, -1.) into +1.) auto sign = +[](double a, double b) { return TMath::Sign(a, b); }; auto signBit = +[](double x) { return std::signbit(x) ? 1.0 : 0.0; }; @@ -188,7 +191,7 @@ std::vector makeTable() F1("TMath::Erfc", +[](double x) { return TMath::Erfc(x); }), F1("TMath::Abs", +[](double x) { return TMath::Abs(x); }, TypeRule::SameAsFirstArg), F1("TMath::Sq", square), - F1("TMath::SignBit", signBit, TypeRule::Int), + F1("TMath::SignBit", signBit, TypeRule::Bool), // two-argument functions F2("pow", +[](double a, double b) { return std::pow(a, b); }), F2("std::pow", +[](double a, double b) { return std::pow(a, b); }), @@ -205,8 +208,8 @@ std::vector makeTable() F2("std::max", stdMax, TypeRule::MinMax), F2("TMath::Max", tmathMax, TypeRule::MinMax), // TFormula shortcut for TMath::Sign - F2("sign", sign, TypeRule::SameAsFirstArg, "TMath::Sign"), - F2("TMath::Sign", sign, TypeRule::SameAsFirstArg), + F2("sign", sign, TypeRule::Sign, "TMath::Sign"), + F2("TMath::Sign", sign, TypeRule::Sign), // zero-argument constants (folded to Op::Const at parse time) F0("TMath::Pi", +[]() { return TMath::Pi(); }), F0("TMath::TwoPi", +[]() { return TMath::TwoPi(); }), diff --git a/roofit/roofitcore/src/RooExprEvaluator.h b/roofit/roofitcore/src/RooExprEvaluator.h index b16ee9ec3c29c..cd1f1d9e424b4 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.h +++ b/roofit/roofitcore/src/RooExprEvaluator.h @@ -28,13 +28,19 @@ namespace RooFormulaFunctions { /// How the C++ result type of a call depends on the argument types. The parser -/// tracks int-ness of subexpressions to reproduce cling's expression typing -/// (in particular to detect integer division, which is not supported). +/// tracks the double/int/bool typing of subexpressions to reproduce cling's +/// expression typing (in particular to detect integer division, which is not +/// supported, and the bool-typed constructs that behave differently). enum class TypeRule : std::uint8_t { Double, ///< result is always double - SameAsFirstArg, ///< result type equals the type of the first argument (abs, sign) - Int, ///< result is an integer type (`int(x)` cast, TMath::SignBit) - MinMax ///< int if both args are int; mixed int/double does not compile in cling + SameAsFirstArg, ///< abs: result type follows the first argument (bool promotes to int) + Int, ///< result is int (the `int(x)` functional cast) + Bool, ///< result is bool (TMath::SignBit) + Sign, ///< sign/TMath::Sign: result type follows the first argument; a bool + ///< first argument is rejected, because cling resolves it to the + ///< generic TMath::Sign template returning bool -- not copysign + MinMax ///< result type is the common argument type; mixed argument types + ///< (int/double or bool/int) do not compile in cling }; struct Entry { diff --git a/roofit/roofitcore/src/RooFormulaParser.cxx b/roofit/roofitcore/src/RooFormulaParser.cxx index 9635bab231320..9ce2d3b650895 100644 --- a/roofit/roofitcore/src/RooFormulaParser.cxx +++ b/roofit/roofitcore/src/RooFormulaParser.cxx @@ -27,13 +27,25 @@ * right-associative exponentiation binding tighter than unary minus, whose * right-hand side may carry one leading sign (see * TFormula::HandleExponentiation). - * - cling's expression typing is tracked (integer literals and bool-valued - * operators are `int` in C++): integer division like `1/2` or `(x>0)/2` - * truncates in cling, so such expressions are not supported here and fall - * back. Similarly min/max with mixed int/double arguments does not compile - * in cling at all and is rejected. + * - cling's expression typing is tracked as double/int/bool: integer + * division like `1/2` or `(x>0)/2` truncates in cling, so such expressions + * are not supported here and fall back. Int-typed constant subexpressions + * are folded in int64 at parse time, and any intermediate leaving the + * int32 range falls back (cling's int arithmetic would wrap around); an + * integer literal too large for int32 (which is long or unsigned in C++, + * not int) falls back as well. min/max with mixed argument types + * (int/double or bool/int) does not compile in cling at all and is + * rejected, and sign()/TMath::Sign with a bool-typed first argument + * resolves to the generic template returning bool (not copysign) and is + * rejected too. * - `%` on doubles does not compile in cling, so TFormula formulas using it * are invalid today; it is not part of this grammar either. + * - Several textual TFormula constructs that are invalid or surprising today + * are kept out of the dialect so that they keep behaving as before (see + * the FallbackTriggers test): the `++` linear-combination separator, runs + * of three or more `-`, bare chained comparisons (cling compiles with + * -Wparentheses as an error), and `^` with a sign on a parenthesized + * exponent (TFormula's rewrite distributes the sign into the group). * - `&&` and `||` do not short-circuit and `?:` evaluates both branches, so * that the scalar and a future vectorized path behave identically. The * selected/combined values are unchanged; this is only observable if @@ -46,6 +58,7 @@ #include #include #include +#include #include #include #include @@ -86,6 +99,7 @@ enum class Tok : std::uint8_t { struct Token { Tok kind = Tok::End; double value = 0.0; ///< Tok::Number + long long intValue = 0; ///< Tok::Number with isInt: the exact integer value bool isInt = false; ///< Tok::Number: literal has integer type in C++ std::uint32_t varIndex = 0; ///< Tok::Var std::string_view text; ///< Tok::Ident / Tok::Number: raw spelling @@ -181,7 +195,13 @@ class Tokenizer { const long long v = std::strtoll(begin, &end, 0); if (errno == ERANGE) return fail("integer literal out of range"); + // An integer literal too large for int32 does not have type int in + // C++ (a decimal one becomes long, a hex/octal one unsigned int), + // which is not the int typing tracked here. Fall back. + if (v > std::numeric_limits::max()) + return fail("integer literal does not fit in int"); tok.value = static_cast(v); + tok.intValue = v; tok.isInt = true; } if (end == begin) @@ -267,8 +287,31 @@ class Tokenizer { Tok kind; std::size_t len = 1; switch (c) { - case '+': kind = Tok::Plus; break; - case '-': kind = Tok::Minus; break; + case '+': + // A textually adjacent `++` is TFormula's linear-combination + // separator (TLinearFitter syntax with one fit parameter per part), + // not an addition. `+ +` with whitespace is ordinary addition. + if (c2 == '+') + return fail("'++' is TFormula's linear-combination separator"); + kind = Tok::Plus; + break; + case '-': { + // TFormula rewrites a double negation `--` (also with whitespace in + // between) so that cling accepts it, but a run of three or more `-` + // survives as a `--` pre-decrement in the generated code, which does + // not compile. Keep such formulas failing (fall back). + std::size_t j = i + 1; + int run = 1; + while (j < n && (_s[j] == '-' || std::isspace(static_cast(_s[j])))) { + if (_s[j] == '-') + ++run; + ++j; + } + if (run >= 3) + return fail("three or more consecutive '-' are invalid in TFormula"); + kind = Tok::Minus; + break; + } case '*': if (c2 == '*') { // TFormula rewrites `**` to `^` kind = Tok::Caret; @@ -384,10 +427,28 @@ class Parser { { } - /// C++ int-vs-double typing of a subexpression, tracked to reject - /// constructs whose cling semantics double arithmetic cannot reproduce. + /// C++ typing of a subexpression (double vs int vs bool), tracked to + /// reject constructs whose cling semantics double arithmetic cannot + /// reproduce: truncating integer division, min/max with mixed argument + /// types, and the bool-typed constructs with non-arithmetic behavior. struct ExprInfo { - bool isInt = false; + enum class Type : std::uint8_t { + Double, + Int, + Bool + }; + Type type = Type::Double; + /// Int-typed constant subexpressions are folded in 64-bit arithmetic at + /// parse time: cling evaluated them in (wrapping) int32 arithmetic, so + /// any int-typed constant intermediate leaving the int32 range makes + /// this evaluator's double arithmetic diverge ("100000*100000" is + /// 1410065408 in cling, not 1e10) and must fall back. Non-constant + /// int-typed intermediates (reachable through an int(x) cast or a + /// promoted bool subexpression) are not tracked. + bool isIntConst = false; + long long intConstValue = 0; + /// Whether the C++ type is an integral type (int or bool). + bool isIntegral() const { return type != Type::Double; } }; bool run(Program &prog) @@ -462,6 +523,11 @@ class Parser { int &_d; }; + static bool fitsInInt32(long long v) + { + return v >= std::numeric_limits::min() && v <= std::numeric_limits::max(); + } + bool fail(std::string msg) { if (_error.empty()) @@ -503,8 +569,16 @@ class Parser { if (!parseTernary(right)) return false; emit(Op::Select); - // The C++ type of `c ? a : b` is int only if both branches are int. - out.isInt = left.isInt && right.isInt; + // The C++ type of `c ? a : b` is the common type of the branches. + if (left.type == ExprInfo::Type::Double || right.type == ExprInfo::Type::Double) { + out.type = ExprInfo::Type::Double; + } else if (left.type == ExprInfo::Type::Bool && right.type == ExprInfo::Type::Bool) { + out.type = ExprInfo::Type::Bool; + } else { + out.type = ExprInfo::Type::Int; + } + // Not a constant (the branch values themselves stay in int32 range). + out.isIntConst = false; return true; } @@ -513,10 +587,19 @@ class Parser { { if (!parseUnary(out)) return false; + // Whether the expression accumulated so far is a bare (unparenthesized) + // relational comparison: cling compiles TFormula code with clang's + // -Wparentheses promoted to an error, so a chained comparison like + // `a < b < c` is invalid in TFormula today (chained equality like + // `a == b == c` is accepted, and parenthesized operands are fine). + bool lhsIsBareRelational = false; while (true) { BinOpInfo const *info = findBinOp(peek().kind); if (!info || info->prec < minPrec) break; + if (info->prec == 4 && lhsIsBareRelational) + return fail("chained comparison is invalid in TFormula"); + lhsIsBareRelational = info->prec == 4; next(); ExprInfo rhs; if (!parseBinary(info->prec + 1, rhs)) @@ -525,24 +608,45 @@ class Parser { case Op::Add: case Op::Sub: case Op::Mul: - out.isInt = out.isInt && rhs.isInt; + // integral operands (bool promotes to int) give an int result + if (out.isIntegral() && rhs.isIntegral()) { + out.type = ExprInfo::Type::Int; + // Fold int-typed constants in 64-bit arithmetic; leaving the + // int32 range means cling's int arithmetic wrapped around, + // which is not reproduced here. Fall back. (Operands are + // within int32 range, so the int64 fold cannot overflow.) + if (out.isIntConst && rhs.isIntConst) { + const long long l = out.intConstValue; + const long long r = rhs.intConstValue; + out.intConstValue = info->op == Op::Add ? l + r : info->op == Op::Sub ? l - r : l * r; + if (!fitsInInt32(out.intConstValue)) + return fail("integer constant expression overflows int in cling"); + } else { + out.isIntConst = false; + } + } else { + out.type = ExprInfo::Type::Double; + out.isIntConst = false; + } emit(info->op); // cling would compute `int * int` in integer arithmetic, where // e.g. (-1) * 0 is +0 and not the -0.0 of double arithmetic. - if (out.isInt && info->op == Op::Mul) + if (out.type == ExprInfo::Type::Int && info->op == Op::Mul) emit(Op::IntNorm); break; case Op::Div: // `1/2` is a truncating integer division in cling. Not supported; // fall back to TFormula so the behavior is unchanged. - if (out.isInt && rhs.isInt) + if (out.isIntegral() && rhs.isIntegral()) return fail("integer division has truncating semantics in TFormula/cling"); - out.isInt = false; + out.type = ExprInfo::Type::Double; + out.isIntConst = false; emit(info->op); break; default: // comparisons and logical operators: C++ result type is bool - out.isInt = true; + out.type = ExprInfo::Type::Bool; + out.isIntConst = false; emit(info->op); break; } @@ -557,21 +661,36 @@ class Parser { if (_depth > kMaxRecursionDepth) return fail("expression too deeply nested"); switch (peek().kind) { - case Tok::Plus: next(); return parseUnary(out); // unary plus: no-op, type preserved + case Tok::Plus: + next(); + if (!parseUnary(out)) + return false; + // unary plus: no-op on the value, but bool promotes to int + if (out.type == ExprInfo::Type::Bool) + out.type = ExprInfo::Type::Int; + return true; case Tok::Minus: next(); if (!parseUnary(out)) return false; emit(Op::Neg); - if (out.isInt) + if (out.isIntegral()) { emit(Op::IntNorm); // cling: -(int)0 is +0, not -0.0 + out.type = ExprInfo::Type::Int; + if (out.isIntConst) { + out.intConstValue = -out.intConstValue; + if (!fitsInInt32(out.intConstValue)) // -(INT_MIN) overflows in cling + return fail("integer constant expression overflows int in cling"); + } + } return true; case Tok::Not: next(); if (!parseUnary(out)) return false; emit(Op::Not); - out.isInt = true; + out.type = ExprInfo::Type::Bool; + out.isIntConst = false; return true; default: return parsePower(out); } @@ -611,7 +730,8 @@ class Parser { } else { emit(Op::Pow); } - out.isInt = false; // pow() and TMath::Sq(Double_t) return double + out.type = ExprInfo::Type::Double; // pow() and TMath::Sq(Double_t) return double + out.isIntConst = false; return true; } @@ -632,6 +752,11 @@ class Parser { haveSign = true; negate = true; } + // TFormula's textual rewrite pushes an explicit sign into a + // parenthesized exponent group, onto only its first term: `x^-(a+b)` + // compiles as pow(x,-(a)+b) today. Do not reproduce that; fall back. + if (haveSign && peek().kind == Tok::LParen) + return fail("'^' with a sign on a parenthesized exponent is broken in TFormula"); const std::size_t operandStart = _pos; if (!parsePower(out)) return false; @@ -641,8 +766,15 @@ class Parser { _tokens[operandStart].text == "2"; if (negate) { emit(Op::Neg); - if (out.isInt) + if (out.isIntegral()) { emit(Op::IntNorm); + out.type = ExprInfo::Type::Int; + if (out.isIntConst) { + out.intConstValue = -out.intConstValue; + if (!fitsInInt32(out.intConstValue)) + return fail("integer constant expression overflows int in cling"); + } + } } return true; } @@ -655,7 +787,9 @@ class Parser { case Tok::Number: { Token const &tok = next(); emit(Op::Const, 0, tok.value); - out.isInt = tok.isInt; + out.type = tok.isInt ? ExprInfo::Type::Int : ExprInfo::Type::Double; + out.isIntConst = tok.isInt; + out.intConstValue = tok.intValue; return true; } case Tok::Var: { @@ -664,7 +798,8 @@ class Parser { return fail("formula references x[" + std::to_string(tok.varIndex) + "] but fewer variables were provided"); _used[tok.varIndex] = true; emit(Op::Var, tok.varIndex); - out.isInt = false; + out.type = ExprInfo::Type::Double; + out.isIntConst = false; return true; } case Tok::LParen: { @@ -717,18 +852,33 @@ class Parser { } using RooFormulaFunctions::TypeRule; + using Type = ExprInfo::Type; switch (entry->rule) { - case TypeRule::Double: out.isInt = false; break; - case TypeRule::SameAsFirstArg: out.isInt = argInfo[0].isInt; break; - case TypeRule::Int: out.isInt = true; break; + case TypeRule::Double: out.type = Type::Double; break; + case TypeRule::SameAsFirstArg: + // abs(bool) resolves to abs(int) in cling: bool promotes to int + out.type = argInfo[0].type == Type::Bool ? Type::Int : argInfo[0].type; + break; + case TypeRule::Int: out.type = Type::Int; break; + case TypeRule::Bool: out.type = Type::Bool; break; + case TypeRule::Sign: + // With a bool first argument, cling resolves TMath::Sign to the + // generic template returning bool: Sign(true, -1.) is +1 there, not + // the -1 of copysign. Fall back rather than reproducing that. + if (argInfo[0].type == Type::Bool) + return fail("'" + name + "' with a bool-typed first argument is not copysign in cling"); + out.type = argInfo[0].type; + break; case TypeRule::MinMax: - // e.g. std::min(x, 3) with double x and int 3 does not compile in - // cling, so such formulas are invalid in TFormula today. Keep it so. - if (argInfo[0].isInt != argInfo[1].isInt) - return fail("'" + name + "' with mixed int/double arguments is invalid in TFormula"); - out.isInt = argInfo[0].isInt; + // e.g. std::min(x, 3) with double x and int 3 (or a bool/int mix) + // does not compile in cling, so such formulas are invalid in + // TFormula today. Keep it so. + if (argInfo[0].type != argInfo[1].type) + return fail("'" + name + "' with mixed argument types is invalid in TFormula"); + out.type = argInfo[0].type; break; } + out.isIntConst = false; // call results are not constant-folded switch (nArgs) { case 0: @@ -740,8 +890,8 @@ class Parser { case 3: emit(Op::Call3, index); break; case 4: emit(Op::Call4, index); break; } - if (out.isInt) - emit(Op::IntNorm); // int-valued calls cannot yield -0.0 in cling + if (out.isIntegral()) + emit(Op::IntNorm); // integer-valued calls cannot yield -0.0 in cling return true; } diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index b0c0d0da706d7..23b3bd663f7f4 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -64,7 +64,9 @@ ROOT_ADD_GTEST(testRooAbsPdf testRooAbsPdf.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooDataSet testRooDataSet.cxx LIBRARIES Tree RooFitCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/dataSet_with_errors_6_26_10.root) ROOT_ADD_GTEST(testRooFormula testRooFormula.cxx LIBRARIES RooFitCore ROOT::TestSupport) -ROOT_ADD_GTEST(testRooFormulaEvaluator testRooFormulaEvaluator.cxx LIBRARIES RooFitCore Hist MathCore) +ROOT_ADD_GTEST(testRooFormulaEvaluator testRooFormulaEvaluator.cxx testRooFormulaEvaluatorIO.cxx + LIBRARIES RooFitCore Hist MathCore + COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testRooFormulaEvaluator_legacy_ws.root) ROOT_ADD_GTEST(testProxiesAndCategories testProxiesAndCategories.cxx LIBRARIES RooFitCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testProxiesAndCategories_1.root diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx index 682d0eb8d9c5c..95b26f0f81cf3 100644 --- a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -372,6 +372,48 @@ TEST(RooFormulaEvaluator, FallbackTriggers) EXPECT_FALSE(astParses("pow(x[0],2^3)")); EXPECT_FALSE(astParses("x[0]>0?1:x[0]^2")); EXPECT_TRUE(astParses("x[0]>0?1:(x[0]^2)")); + // `^` with an explicit sign on a parenthesized exponent: TFormula's + // textual rewrite pushes the sign onto only the first term inside the + // group (`x^-(a+b)` compiles as pow(x,-(a)+b) today), so fall back + EXPECT_FALSE(astParses("x[0]^-(x[0]+1)")); + EXPECT_FALSE(astParses("x[0]^-(2)")); + EXPECT_FALSE(astParses("x[0]^+(x[0]*2)")); + EXPECT_TRUE(astParses("x[0]^(-x[0]-1)")); // sign inside the group is fine + EXPECT_TRUE(astParses("x[0]^-2.5")); + EXPECT_TRUE(astParses("x[0]^-sin(x[0])")); + // cling resolves TMath::Sign with a bool first argument to the generic + // template, which returns bool -- not copysign + EXPECT_FALSE(astParses("sign(x[0]>1, x[0])")); + EXPECT_FALSE(astParses("TMath::Sign(!x[0], x[0])")); + EXPECT_FALSE(astParses("sign(TMath::SignBit(x[0]), x[0])")); + EXPECT_TRUE(astParses("sign(1, x[0])")); // an int first argument is copysign-like + // bool/int mixes in min/max do not compile in cling either + EXPECT_FALSE(astParses("min(x[0]>1, 2)")); + EXPECT_TRUE(astParses("min(x[0]>1, x[0]>2)")); + EXPECT_TRUE(astParses("abs(x[0]>1)")); // abs(bool) promotes to int and works + // cling compiles the TFormula code with -Wparentheses promoted to an + // error, so a bare chained comparison is invalid in TFormula today + EXPECT_FALSE(astParses("x[0]1")); + EXPECT_FALSE(astParses("1+x[0]<2<3")); + EXPECT_TRUE(astParses("(x[0]1")); + // a textually adjacent `++` is TFormula's linear-combination separator + // (one fit parameter per part), not an addition + EXPECT_FALSE(astParses("2++3")); + EXPECT_FALSE(astParses("x[0]++x[1]")); + EXPECT_FALSE(astParses("++x[0]")); + EXPECT_TRUE(astParses("2+ +3")); // with whitespace it is an ordinary addition + // runs of three or more `-` (with or without whitespace) survive TFormula's + // double-negation rewrite as a `--` pre-decrement, which cling rejects + EXPECT_FALSE(astParses("2---3")); + EXPECT_FALSE(astParses("---x[0]")); + EXPECT_FALSE(astParses("x[0]- - -x[1]", 2)); + EXPECT_TRUE(astParses("2--3")); + EXPECT_TRUE(astParses("x[0]- -x[1]", 2)); + EXPECT_TRUE(astParses("-(-(-x[0]))")); // TFormula parameters and parametrized shortcuts EXPECT_FALSE(astParses("[0]+x[0]")); EXPECT_FALSE(astParses("x[0]+[0]*0xaf")); @@ -405,6 +447,63 @@ TEST(RooFormulaEvaluator, FallbackTriggers) EXPECT_TRUE(astParses(shallow)); } +// Int-typed constant arithmetic that leaves the int32 range wrapped around in +// cling's int arithmetic ("100000*100000" is 1410065408 there, not 1e10). +// That is not reproduced in double arithmetic: such formulas must fall back +// to the TFormula backend, so their values are unchanged. Integer literals +// too large for int32 have type long (decimal) or unsigned int (hex/octal) in +// C++ -- not the int typing tracked by the parser -- and fall back too. +TEST(RooFormulaEvaluator, IntOverflowFallback) +{ + EXPECT_FALSE(astParses("100000*100000")); + EXPECT_FALSE(astParses("x[0]*(100000*100000)")); + EXPECT_FALSE(astParses("2000000000+2000000000")); + EXPECT_FALSE(astParses("2147483647+1")); + EXPECT_FALSE(astParses("0-2000000000-2000000000")); + EXPECT_FALSE(astParses("-(0-2147483647-1)")); // negating INT_MIN overflows + EXPECT_TRUE(astParses("100000*100000.")); // double arithmetic is fine + EXPECT_TRUE(astParses("100000.*100000")); + EXPECT_EQ(astVal("2147483647*1"), 2147483647.0); + EXPECT_EQ(astVal("0-2147483647-1"), -2147483648.0); // INT_MIN itself is in range + // literals out of int32 range + EXPECT_FALSE(astParses("3000000000")); + EXPECT_FALSE(astParses("0x80000000")); + EXPECT_FALSE(astParses("min(3000000000,2)")); + EXPECT_TRUE(astParses("2147483647")); + EXPECT_TRUE(astParses("0x7fffffff")); + + ScopedBackendEnv env{nullptr}; + RooRealVar x("x", "x", 5.0); + + // The fallback must reproduce cling's wrapped value exactly. + RooArgList vars{x}; + auto f = RooFormulaUtils::makeFormulaEvaluator("f", "x*(100000*100000)", vars); + EXPECT_FALSE(isAstBackend(*f)); + TFormula ref("ref", "x*(100000*100000)", /*addToGlobList=*/false); + ASSERT_TRUE(ref.IsValid()); + double xv = 5.0; + EXPECT_TRUE(sameBits(RooFormulaUtils::evalFormula(*f, vars), ref.EvalPar(&xv))); + + // A long-typed literal on its own is valid in cling with the same value; + // the fallback keeps such formulas working. + auto g = RooFormulaUtils::makeFormulaEvaluator("g", "3000000000+0*x", vars); + EXPECT_FALSE(isAstBackend(*g)); + EXPECT_EQ(RooFormulaUtils::evalFormula(*g, vars), 3000000000.0); + + // min(long, int) does not compile in cling, so this formula threw at + // construction before the JIT-free backend existed. It must still throw + // instead of silently yielding 2. + { + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TFormula::InputFormulaIntoCling", "Error compiling formula expression in Cling", + false); + diags.requiredDiag(kError, "TFormula::ProcessFormula", " is invalid", false); + diags.optionalDiag(kError, "prepareMethod", "Can't compile function TFormula", false); + diags.optionalDiag(kError, "cling", "no matching function", false); + EXPECT_THROW(RooFormulaUtils::makeFormulaEvaluator("h", "min(3000000000,2)+x", vars), std::runtime_error); + } +} + // "0x1e+2" is one single (invalid) pp-number in C++, not 0x1e + 2: cling // refused it, so formula construction threw. The lexer must not split // it into two tokens; such formulas keep failing via the TFormula fallback. @@ -903,6 +1002,11 @@ TEST(RooFormulaEvaluator, DifferentialCorpus) for (auto const &f : fallbacks) { std::cout << " fallback: " << f << "\n"; } + RecordProperty("CorpusSize", nTotal); + RecordProperty("CorpusOnAstPath", nAst); + RecordProperty("CorpusCoveragePercent", static_cast(std::round(100. * coverage))); + // Hard floor only: the detailed number is reported above. A drop below 90% + // means the allow-list lost something that real-world formulas need. EXPECT_GE(coverage, 0.9); } @@ -1033,3 +1137,395 @@ TEST(RooFormulaEvaluator, EmitCppLocaleIndependent) EXPECT_NE(expr.find("1.25"), std::string::npos) << expr; EXPECT_EQ(expr.find(','), std::string::npos) << expr; } + +namespace { + +/// A randomly generated expression in the processed `x[i]` dialect, together +/// with the C++ type (double, int, or bool) its cling compilation would have. +struct RandomExpr { + enum class Type : std::uint8_t { + Double, + Int, + Bool + }; + std::string text; + Type type = Type::Double; + bool isIntegral() const { return type != Type::Double; } +}; + +/// Generates random well-formed expressions over exactly the grammar the +/// JIT-free parser supports: all binary and unary operators including `^` and +/// `**`, comparisons, logical operators, the ternary operator, and every +/// function spelling in the RooFormulaFunctions allow-list (sampled directly +/// from the table, so new entries are covered automatically, including the +/// zero-argument constants and the multi-arity TMath::Gaus). +/// +/// The generator composes strings level by level along the C++ operator +/// precedence, so the string parses back to the generated structure and the +/// tracked double/int/bool typing is that of the actual parse. The typing is +/// used to steer around the deliberately unsupported constructs (truncating +/// integer division, min/max with mixed argument types, sign() with a +/// bool-typed first argument) by inserting a `1.0 *` factor. The textual +/// pitfalls of the dialect are avoided structurally: `^` expressions are +/// always parenthesized (a `^` operand adjacent to `,` or `:` is invalid in +/// TFormula), a signed exponent never starts with `(` (TFormula's rewrite +/// distributes the sign into the group), comparisons are not chained, and +/// stacked signs are parenthesized (`++` is TFormula's linear-combination +/// separator and runs of three `-` are invalid). Everything the generator +/// produces must therefore parse on the AST path *and* be a valid TFormula, +/// so the differential test can require bitwise agreement with no skips. +class RandomExprGenerator { +public: + RandomExprGenerator(unsigned int seed, unsigned int nVars) : _rng{seed}, _nVars{nVars} {} + + RandomExpr gen(int depth) { return genTernary(depth); } + +private: + double chance() { return std::uniform_real_distribution{0., 1.}(_rng); } + int pick(int n) { return std::uniform_int_distribution{0, n - 1}(_rng); } + + /// Turn an int- or bool-typed operand into a double-typed one with the + /// same value. + RandomExpr forceDouble(RandomExpr e) + { + if (e.isIntegral()) { + e.text = "(1.0 * (" + e.text + "))"; + e.type = RandomExpr::Type::Double; + } + return e; + } + + RandomExpr genTernary(int depth) + { + if (depth > 0 && chance() < 0.10) { + RandomExpr c = genBinary(1, depth - 1); + RandomExpr a = genTernary(depth - 1); + RandomExpr b = genTernary(depth - 1); + RandomExpr out; + out.text = c.text + " ? " + a.text + " : " + b.text; + if (a.type == RandomExpr::Type::Double || b.type == RandomExpr::Type::Double) { + out.type = RandomExpr::Type::Double; + } else if (a.type == RandomExpr::Type::Bool && b.type == RandomExpr::Type::Bool) { + out.type = RandomExpr::Type::Bool; + } else { + out.type = RandomExpr::Type::Int; + } + return out; + } + return genBinary(1, depth); + } + + /// Precedence levels: 1 `||`, 2 `&&`, 3 `== !=`, 4 `< <= > >=`, 5 `+ -`, + /// 6 `* /`; operands of a level come from the next-tighter level, so no + /// parentheses are needed to reproduce the intended structure. + RandomExpr genBinary(int level, int depth) + { + if (level > 6) + return genUnary(depth); + RandomExpr lhs = genBinary(level + 1, depth); + static constexpr double probs[7] = {0.0, 0.06, 0.06, 0.08, 0.10, 0.35, 0.35}; + while (depth > 0 && chance() < probs[level]) { + --depth; + RandomExpr rhs = genBinary(level + 1, depth); + const char *op = nullptr; + switch (level) { + case 1: + op = "||"; + lhs.type = RandomExpr::Type::Bool; + break; + case 2: + op = "&&"; + lhs.type = RandomExpr::Type::Bool; + break; + case 3: + op = pick(2) ? "==" : "!="; + lhs.type = RandomExpr::Type::Bool; + break; + case 4: + switch (pick(4)) { + case 0: op = "<"; break; + case 1: op = "<="; break; + case 2: op = ">"; break; + default: op = ">="; break; + } + lhs.type = RandomExpr::Type::Bool; + break; + case 5: + op = pick(2) ? "+" : "-"; + lhs.type = lhs.isIntegral() && rhs.isIntegral() ? RandomExpr::Type::Int : RandomExpr::Type::Double; + break; + case 6: + if (pick(3) == 0) { + // avoid truncating integer division (unsupported: falls back) + if (lhs.isIntegral()) + rhs = forceDouble(rhs); + op = "/"; + lhs.type = RandomExpr::Type::Double; + } else { + op = "*"; + lhs.type = lhs.isIntegral() && rhs.isIntegral() ? RandomExpr::Type::Int : RandomExpr::Type::Double; + } + break; + } + lhs.text += std::string{" "} + op + " " + rhs.text; + // no bare chained comparison (`a < b < c`): invalid in TFormula, + // where cling compiles with -Wparentheses promoted to an error + if (level == 4) + break; + } + return lhs; + } + + RandomExpr genUnary(int depth) + { + if (depth > 0 && chance() < 0.15) { + RandomExpr e = genUnary(depth - 1); + const int which = pick(3); // -, +, ! + if (which == 2) { + e.text = "!(" + e.text + ")"; + e.type = RandomExpr::Type::Bool; + } else { + // parenthesize an operand that starts with a sign or `!`, to + // avoid the `++` / `---` pitfalls (see FallbackTriggers) + if (e.text[0] == '-' || e.text[0] == '+' || e.text[0] == '!') + e.text = "(" + e.text + ")"; + e.text = (which == 0 ? "-" : "+") + e.text; + if (e.isIntegral()) + e.type = RandomExpr::Type::Int; // bool promotes to int + } + return e; + } + return genPower(depth); + } + + /// `^`/`**` exponentiation, always parenthesized as a whole; base and + /// exponent are primaries (with one optional sign on the exponent), + /// matching how the operator appears in real formulas and staying within + /// what TFormula's textual pow() rewrite scans correctly. + RandomExpr genPower(int depth) + { + if (depth > 0 && chance() < 0.12) { + RandomExpr base = genPrimary(depth - 1); + std::string sign; + if (chance() < 0.3) + sign = pick(2) ? "-" : "+"; + RandomExpr exponent = genPrimary(depth - 1); + // a signed exponent must not start with `(`: TFormula's textual + // rewrite distributes the sign into the group (falls back) + if (!sign.empty() && exponent.text[0] == '(') + sign.clear(); + const char *op = pick(4) == 0 ? "**" : "^"; + RandomExpr out; + out.text = "(" + base.text + op + sign + exponent.text + ")"; + out.type = RandomExpr::Type::Double; + return out; + } + return genPrimary(depth); + } + + RandomExpr genPrimary(int depth) + { + const double r = chance(); + if (depth <= 0 || r < 0.4) { + if (chance() < 0.55) { + RandomExpr out; + out.text = "x[" + std::to_string(pick(_nVars)) + "]"; + return out; + } + return genLiteral(); + } + if (r < 0.6) { + RandomExpr e = genTernary(depth - 1); + e.text = "(" + e.text + ")"; + return e; + } + return genCall(depth - 1); + } + + RandomExpr genLiteral() + { + static const char *const kIntLiterals[] = {"0", "1", "2", "3", "7", "42", "0x1f"}; + static const char *const kDoubleLiterals[] = {"0.5", "1.5", "2.5", "3.360779", "0.25", ".5", + "1.", "1e-3", "2e+2", "1e300", "1e-300", "0.1", + "1e30", "13.7", "1.5e-8", "2.0", "6.62607015e-34"}; + RandomExpr out; + if (chance() < 0.4) { + out.text = kIntLiterals[pick(std::end(kIntLiterals) - std::begin(kIntLiterals))]; + out.type = RandomExpr::Type::Int; + } else { + out.text = kDoubleLiterals[pick(std::end(kDoubleLiterals) - std::begin(kDoubleLiterals))]; + } + return out; + } + + /// A call to a random entry of the actual allow-list table. + RandomExpr genCall(int depth) + { + auto const *tab = RooFormulaFunctions::table(); + auto const &entry = tab[pick(static_cast(RooFormulaFunctions::tableSize()))]; + + RandomExpr args[4]; + for (unsigned int i = 0; i < entry.arity; ++i) { + args[i] = genTernary(depth); + } + + using RooFormulaFunctions::TypeRule; + using Type = RandomExpr::Type; + Type type = Type::Double; + switch (entry.rule) { + case TypeRule::Double: break; + case TypeRule::SameAsFirstArg: + if (entry.arity >= 1) + type = args[0].type == Type::Bool ? Type::Int : args[0].type; + break; + case TypeRule::Int: type = Type::Int; break; + case TypeRule::Bool: type = Type::Bool; break; + case TypeRule::Sign: + // sign() with a bool first argument is not copysign in cling (falls back) + if (args[0].type == Type::Bool) + args[0] = forceDouble(args[0]); + type = args[0].type; + break; + case TypeRule::MinMax: + // mixed argument types do not compile in cling (falls back) + if (args[0].type != args[1].type) { + args[0] = forceDouble(args[0]); + args[1] = forceDouble(args[1]); + } + type = args[0].type; + break; + } + + RandomExpr out; + out.type = type; + out.text = std::string{entry.name} + "("; + for (unsigned int i = 0; i < entry.arity; ++i) { + if (i > 0) + out.text += ", "; + out.text += args[i].text; + } + out.text += ")"; + return out; + } + + std::mt19937 _rng; + unsigned int _nVars = 0; +}; + +constexpr unsigned int kRandomExprSeed = 20260827u; // fixed: failures must reproduce +constexpr unsigned int kRandomExprVars = 3; + +/// Fill the input vector for the given trial: the first trials use edge +/// values (0, +-1, very large, very small, negatives that drive sqrt/log/pow +/// into NaN territory), later ones random values from the given generator. +void fillRandomInputs(int trial, double *pars, std::mt19937 &rng) +{ + switch (trial) { + case 0: + for (unsigned int j = 0; j < kRandomExprVars; ++j) + pars[j] = 0.0; + break; + case 1: + for (unsigned int j = 0; j < kRandomExprVars; ++j) + pars[j] = j == 1 ? -1.0 : 1.0; + break; + case 2: + pars[0] = 1e300; + pars[1] = 1e-300; + pars[2] = -1.0; + break; + case 3: + pars[0] = -2.5; + pars[1] = -1e300; + pars[2] = -1e-300; + break; + default: { + const double scale = trial % 2 ? 3.0 : 50.0; + std::uniform_real_distribution dist{-scale, scale}; + for (unsigned int j = 0; j < kRandomExprVars; ++j) + pars[j] = dist(rng); + break; + } + } +} + +} // namespace + +// Random-expression differential campaign: several hundred generated +// expressions, each evaluated through both backends on several input vectors +// including edge cases. Results must agree bitwise (NaN counts as equal to +// NaN); there is no tolerance to widen. The seed is fixed, so a failure +// reproduces exactly; to investigate one, print `expr.text` and the +// hexfloat values from the failure message. +TEST(RooFormulaEvaluator, RandomExpressionDifferential) +{ + constexpr int nExprs = 500; + constexpr int nTrials = 8; + + RandomExprGenerator gen{kRandomExprSeed, kRandomExprVars}; + std::mt19937 inputRng{987654u}; + + for (int iExpr = 0; iExpr < nExprs; ++iExpr) { + const RandomExpr expr = gen.gen(4); + + std::string error; + auto prog = RooFormulaParser::compile(expr.text, kRandomExprVars, &error); + ASSERT_TRUE(prog) << "generated expression unexpectedly failed to parse: " << expr.text << "\n error: " << error; + RooExprEvaluator ast{prog}; + + TFormula ref("ref", expr.text.c_str(), /*addToGlobList=*/false); + ASSERT_TRUE(ref.IsValid()) << "generated expression is invalid in TFormula (the generator must only produce " + "expressions valid in both dialects): " + << expr.text; + + for (int trial = 0; trial < nTrials; ++trial) { + double pars[kRandomExprVars]; + fillRandomInputs(trial, pars, inputRng); + const double a = ast.eval(pars); + const double t = ref.EvalPar(pars); + EXPECT_TRUE(sameBits(a, t)) << "AST and TFormula backends disagree for: " << expr.text + << "\n inputs: " << pars[0] << " " << pars[1] << " " << pars[2] + << "\n ast = " << std::hexfloat << a << " tformula = " << t << std::defaultfloat; + } + } +} + +// Emitted-code agreement on random expressions (the corpus counterpart is +// EmittedCppDifferential): emit the C++ for a sample of the same generated +// expression stream, compile it with the interpreter, and require bitwise +// agreement with AST evaluation. Uses the same seed as +// RandomExpressionDifferential, so this covers a prefix of the same +// expressions. +TEST(RooFormulaEvaluator, EmittedCppRandomExpressions) +{ + constexpr int nExprs = 150; + constexpr int nTrials = 6; + + auto varName = [](unsigned int i) { return "v[" + std::to_string(i) + "]"; }; + + RandomExprGenerator gen{kRandomExprSeed, kRandomExprVars}; + std::mt19937 inputRng{192837u}; + + for (int iExpr = 0; iExpr < nExprs; ++iExpr) { + const RandomExpr expr = gen.gen(4); + + auto prog = RooFormulaParser::compile(expr.text, kRandomExprVars); + ASSERT_TRUE(prog) << expr.text; + RooExprEvaluator ast{prog}; + + const std::string emitted = ast.emitCpp(varName); + ASSERT_FALSE(emitted.empty()) << expr.text; + EmittedFunc fn = compileEmitted(emitted); + ASSERT_NE(fn, nullptr) << "emitted C++ failed to compile for: " << expr.text << "\n emitted: " << emitted; + + for (int trial = 0; trial < nTrials; ++trial) { + double pars[kRandomExprVars]; + fillRandomInputs(trial, pars, inputRng); + const double a = ast.eval(pars); + const double c = fn(pars); + EXPECT_TRUE(sameBits(a, c)) << "AST and emitted C++ disagree for: " << expr.text << "\n emitted: " << emitted + << "\n inputs: " << pars[0] << " " << pars[1] << " " << pars[2] + << "\n ast = " << std::hexfloat << a << " emitted = " << c << std::defaultfloat; + } + } +} diff --git a/roofit/roofitcore/test/testRooFormulaEvaluatorIO.cxx b/roofit/roofitcore/test/testRooFormulaEvaluatorIO.cxx new file mode 100644 index 0000000000000..1090846267cce --- /dev/null +++ b/roofit/roofitcore/test/testRooFormulaEvaluatorIO.cxx @@ -0,0 +1,296 @@ +// Workspace-I/O tests for the JIT-free RooFormula evaluation backend: +// round-tripping formulas written through the AST path, reading a legacy +// workspace written by a pre-AST (TFormula-backed) build, and checking that +// both backends persist an identical on-disk representation. +// Author: Jonas Rembser, CERN 2026 + +#include "../src/RooFormulaUtils.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace { + +/// Bitwise comparison; any two NaNs count as equal. +bool sameBits(double a, double b) +{ + if (std::isnan(a) && std::isnan(b)) + return true; + return std::memcmp(&a, &b, sizeof(double)) == 0; +} + +/// Set or unset ROOFIT_FORMULA_BACKEND for the lifetime of the object, +/// resetting the backend's read-once cache on both ends. +class ScopedBackendEnv { +public: + ScopedBackendEnv(const char *value) + { + if (const char *old = std::getenv("ROOFIT_FORMULA_BACKEND")) + _old = old; + if (value) + setenv("ROOFIT_FORMULA_BACKEND", value, /*overwrite=*/1); + else + unsetenv("ROOFIT_FORMULA_BACKEND"); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + ~ScopedBackendEnv() + { + if (_old.empty()) + unsetenv("ROOFIT_FORMULA_BACKEND"); + else + setenv("ROOFIT_FORMULA_BACKEND", _old.c_str(), 1); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + +private: + std::string _old; +}; + +const char *const kObjectNames[] = {"f1", "f2", "f3", "f4", "f5", "f6", "f7", "p1"}; + +/// Fill the workspace with the same formula content that the legacy fixture +/// testRooFormulaEvaluator_legacy_ws.root was generated from (by a pre-AST, +/// Phase-1-only build whose persisted output is identical to unpatched ROOT): +/// plain named variables, @-references, TMath:: functions, `^`, a ternary, +/// and a category-state formula. +void fillWorkspace(RooWorkspace &ws) +{ + RooRealVar x("x", "x", 0.5, -10, 10); + RooRealVar a("a", "a", 1.1, 0.1, 5); + RooRealVar b("b", "b", 2.3, 0.1, 5); + + RooFormulaVar f1("f1", "simple", "a*x + b", {a, x, b}); + RooFormulaVar f2("f2", "atref", "@0*@1 - @2", {a, x, b}); + RooFormulaVar f3("f3", "funcs", "sqrt(abs(x)) + TMath::Erf(a) + exp(-b)", {x, a, b}); + RooFormulaVar f4("f4", "pow", "x^2 + a^-0.5 + pow(b, 3)", {x, a, b}); + RooFormulaVar f5("f5", "ternary", "x > 0 ? log(1+x) : -x", {x}); + RooFormulaVar f6("f6", "histfactory-shape", "1 + 0.02*a", {a}); + RooGenericPdf p1("p1", "genpdf", "exp(-0.5*((x-a)/b)^2)", {x, a, b}); + + RooCategory c("c", "c"); + c.defineType("sig", 0); + c.defineType("bkg", 1); + RooFormulaVar f7("f7", "catstate", "c == c::sig ? a : b", {c, a, b}); + + ws.import(f1); + ws.import(f2); + ws.import(f3); + ws.import(f4); + ws.import(f5); + ws.import(f6); + ws.import(p1); + ws.import(f7); +} + +/// The persisted formula string (_formExpr) of a RooFormulaVar or +/// RooGenericPdf, read without constructing the evaluation engine. +std::string persistedExpression(RooAbsArg *arg) +{ + if (auto *formulaVar = dynamic_cast(arg)) { + return formulaVar->expression(); + } + if (auto *genericPdf = dynamic_cast(arg)) { + return genericPdf->expression(); + } + return {}; +} + +/// Map of class name -> streamed class version for all TStreamerInfos in the +/// given file. +std::map streamerInfoVersions(TFile &file) +{ + std::map out; + std::unique_ptr infos{file.GetStreamerInfoList()}; + for (TObject *obj : *infos) { + if (auto *info = dynamic_cast(obj)) { + out[info->GetName()] = info->GetClassVersion(); + } + } + return out; +} + +} // namespace + +// Round-trip regression test for the formulaString() persistence landmine: +// RooFormulaVar/RooGenericPdf must persist a non-empty processed formula +// string when the AST backend is active (no TFormula exists to read a title +// from), and evaluate bitwise identically after reading back. +TEST(RooFormulaEvaluatorIO, WorkspaceRoundTrip) +{ + // Force the AST path: with a silent fallback this test could pass without + // testing anything, and a formula that stops parsing shows up as a throw. + ScopedBackendEnv env{"ast"}; + + const char *fileName = "testRooFormulaEvaluatorIO_roundtrip.root"; + + std::map refValues; + std::map refExprs; + { + RooWorkspace ws{"w"}; + fillWorkspace(ws); + for (const char *name : kObjectNames) { + auto *arg = dynamic_cast(ws.arg(name)); + ASSERT_NE(arg, nullptr) << name; + refValues[name] = arg->getVal(); + // after construction, _formExpr holds the processed x[i]-dialect + // string; it must not be empty + refExprs[name] = persistedExpression(arg); + ASSERT_FALSE(refExprs[name].empty()) << name; + } + ASSERT_TRUE(ws.writeToFile(fileName)); // returns true on success + } + + { + TFile file(fileName, "READ"); + ASSERT_FALSE(file.IsZombie()); + auto *ws = file.Get("w"); + ASSERT_NE(ws, nullptr); + for (const char *name : kObjectNames) { + auto *arg = dynamic_cast(ws->arg(name)); + ASSERT_NE(arg, nullptr) << name; + EXPECT_EQ(persistedExpression(arg), refExprs[name]) << name; + EXPECT_TRUE(sameBits(arg->getVal(), refValues[name])) + << name << ": pre-write = " << std::hexfloat << refValues[name] << " read-back = " << arg->getVal() + << std::defaultfloat; + } + } + + gSystem->Unlink(fileName); +} + +namespace { + +// Reference values recorded (with 17 significant digits, so the double +// literals below reproduce them bitwise) by the pre-AST build that wrote +// testRooFormulaEvaluator_legacy_ws.root. That build evaluated through +// TFormula/cling, so these are the historical values that reading the file +// must reproduce exactly. +const std::pair kLegacyReference[] = { + {"f1", 2.8499999999999996}, // a*x + b + {"f2", -1.7499999999999998}, // @0*@1 - @2 + {"f3", 1.687570694483433}, // sqrt(abs(x)) + TMath::Erf(a) + exp(-b) + {"f4", 13.370462589245591}, // x^2 + a^-0.5 + pow(b, 3) + {"f5", 0.40546510810816438}, // x > 0 ? log(1+x) : -x + {"f6", 1.022}, // 1 + 0.02*a + {"f7", 1.1000000000000001}, // c == c::sig ? a : b + {"p1", 0.96654592463371813}, // exp(-0.5*((x-a)/b)^2) +}; + +} // namespace + +// Reading a workspace written by a pre-AST build must reproduce the recorded +// values bitwise, on the AST path (which all persisted strings must reach: +// they are stored in the processed x[i] dialect the parser targets) as well +// as on the TFormula fallback path. +TEST(RooFormulaEvaluatorIO, LegacyWorkspaceValues) +{ + for (const char *backend : {"ast", "tformula"}) { + ScopedBackendEnv env{backend}; + TFile file("testRooFormulaEvaluator_legacy_ws.root", "READ"); + ASSERT_FALSE(file.IsZombie()); + auto *ws = file.Get("w"); + ASSERT_NE(ws, nullptr); + for (auto const &ref : kLegacyReference) { + auto *arg = dynamic_cast(ws->arg(ref.first)); + ASSERT_NE(arg, nullptr) << ref.first; + const double val = arg->getVal(); + EXPECT_TRUE(sameBits(val, ref.second)) + << ref.first << " (backend " << backend << "): expected = " << std::hexfloat << ref.second + << " actual = " << val << std::defaultfloat; + } + } +} + +// Fidelity of the persisted form: writing the same workspace content through +// the AST backend and through the TFormula backend must produce identical +// persisted formula strings and identical streamed class versions -- i.e. the +// on-disk representation does not depend on the evaluation backend, and it +// matches what the pre-AST build wrote. (A true cross-version read test with +// an unpatched ROOT release is complementary to this in-process check.) +TEST(RooFormulaEvaluatorIO, BackendWriteFidelity) +{ + const char *fileNameAst = "testRooFormulaEvaluatorIO_fidelity_ast.root"; + const char *fileNameTF = "testRooFormulaEvaluatorIO_fidelity_tformula.root"; + + { + ScopedBackendEnv env{"ast"}; + RooWorkspace ws{"w"}; + fillWorkspace(ws); + ASSERT_TRUE(ws.writeToFile(fileNameAst)); + } + { + ScopedBackendEnv env{"tformula"}; + RooWorkspace ws{"w"}; + fillWorkspace(ws); + ASSERT_TRUE(ws.writeToFile(fileNameTF)); + } + + TFile fileAst(fileNameAst, "READ"); + TFile fileTF(fileNameTF, "READ"); + TFile fileLegacy("testRooFormulaEvaluator_legacy_ws.root", "READ"); + ASSERT_FALSE(fileAst.IsZombie()); + ASSERT_FALSE(fileTF.IsZombie()); + ASSERT_FALSE(fileLegacy.IsZombie()); + + auto *wsAst = fileAst.Get("w"); + auto *wsTF = fileTF.Get("w"); + auto *wsLegacy = fileLegacy.Get("w"); + ASSERT_NE(wsAst, nullptr); + ASSERT_NE(wsTF, nullptr); + ASSERT_NE(wsLegacy, nullptr); + + for (const char *name : kObjectNames) { + auto *argAst = dynamic_cast(wsAst->arg(name)); + auto *argTF = dynamic_cast(wsTF->arg(name)); + auto *argLegacy = dynamic_cast(wsLegacy->arg(name)); + ASSERT_NE(argAst, nullptr) << name; + ASSERT_NE(argTF, nullptr) << name; + ASSERT_NE(argLegacy, nullptr) << name; + + const std::string exprAst = persistedExpression(argAst); + EXPECT_FALSE(exprAst.empty()) << name; + EXPECT_EQ(exprAst, persistedExpression(argTF)) << name; + EXPECT_EQ(exprAst, persistedExpression(argLegacy)) << name; + + EXPECT_TRUE(sameBits(argAst->getVal(), argTF->getVal())) << name; + EXPECT_TRUE(sameBits(argAst->getVal(), argLegacy->getVal())) << name; + } + + // The streamed class versions must be independent of the backend (no + // schema change of any kind)... + const auto versionsAst = streamerInfoVersions(fileAst); + const auto versionsTF = streamerInfoVersions(fileTF); + const auto versionsLegacy = streamerInfoVersions(fileLegacy); + EXPECT_FALSE(versionsAst.empty()); + EXPECT_EQ(versionsAst, versionsTF); + // ... and every class streamed by the current code must keep the version + // it had in the legacy file. (The legacy file additionally contains a + // pair dictionary entry from the "origName" string + // attributes that the RooFormula-based code stamped on imported formula + // dependents; the current code doesn't write these noise attributes.) + for (auto const &entry : versionsAst) { + auto it = versionsLegacy.find(entry.first); + ASSERT_NE(it, versionsLegacy.end()) << entry.first; + EXPECT_EQ(it->second, entry.second) << entry.first; + } + + gSystem->Unlink(fileNameAst); + gSystem->Unlink(fileNameTF); +} diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator_legacy_ws.root b/roofit/roofitcore/test/testRooFormulaEvaluator_legacy_ws.root new file mode 100644 index 0000000000000000000000000000000000000000..9b17c3d055d70209f3a7d5f78d771f24a4ff4812 GIT binary patch literal 9243 zcmb`tWl$YW)Gd5)C%C&4EVu`EcXtkOu!Flx(BKd>I0U!g8r&_oySqD=Cw1%le%^bl z-mcm`UDMNR_L`}gy{q?f0)t%ufF~dT05Asts0aZ7%A$|4>xV&nSOfG&_b~}y3;+o>=6GraN z%>S(e{0~9_VE*&=!;C&0J^!$V5Buje06@|Be|!`f{}@aDW32HXV+hUv)+4GbJAuuB z&dxI8i2vW{3zY>cD>pYgmj#mv8yhDRE2{-7lPMQ3E7KP?E)EkOvoEGVpoK61APRu^ z=mkFD{(Cw9+W#vG007(YpKG{i|EMwibM5j2${bl04AufW**H6xm;oW&A>4K4;e%iS zXt)3mh%s#=TQwk-(E90|Yat;aA+oqkV?=q6h?clZQ0vwF1Kf}@2^rMc6S!|`u0i^; z1l;kE1B>Osx`}{U1jaLS2{ASeKZli1g2Plm3%xKE5AJk6JAGevixGK9@oijpA9o*b z-k&(!w;zroidRLh>2Tlc39~(8rX!qBi{#HZt zZ3Duz`B|XZOGx@8Q$+oDFlh67JbUE(Kn;+uL;k;%Cd`1mFgnxfr?62{dP$n z(*R}~SZEY7#xxr^%2Y4T&*sSb%}%6rn%e73?^Jy8P1ih_!ODCZpHeNr)itGg*4rXt z>o31%hb@E}M?7EuA_n?TgUmr|22PVXTVLqmnA01pgI;33{XMFKGx>@oBCv}Ro zqIi;=2i?4B($Dr{BoRc^sM<29LY?g@+7Bf?i@>h27&MGNT_27dl@rBaX`eo#a-3!q z+%nJ|%Dq(s&o?wy;~{_(>jcj=d@LYJ@U^igi~6WK81UT>1Z#9&R0(DHeJ(f|;6wcKQPmfQtJ9^D@2(?1bI1GCieX zL>tdy%%N--vlPdAY@6yU(!@}~5uVuG0qL;72JSB32-5T%j}<=4)MOrOwEh(G)V zd~+jhDWf#i_%qWF*lQ%0bvADT^HbRH9fB+uT*=|6$W?nNq{^tTft`qZc#%2ntjM!> zR#3mL^|R|^-~g_Re9GKjC~wTD?dUUhWAV=fO$3n3p#usfEw>z6bB4ag8KrQxUoLV% zkeLp>0-O-+?hmNV-t$hF(93~@>m!AIDj^;m0Azd8CB?#Xfypcb{M`>SJOHWx+$4pPyJG)z8hl4mE-{I$62Xw@ee{#*2Ez=?d$n@y zmQ9D7?^x*c1+#VFuqWdD!B~6jPiK@W3TNSJ;*|OTU4lvsqlIrJ&blhgUn=q53H$-M zwPDOWQ1xy48|Yrz%`2G7yRN;>&(k!U2Gif^=m^`HsbVtsEJZ zcE`54CnR3Kh^IzwF~9YH2Y1i(10}MXmpG5`a03x;)6)1zF9-320&mbds^R0Aq-QT}^F0 z$!tOPHbC=_gR{#=C>y{BFlc;)vQzy}(L`FF?`~&4caJAL4Q90Cq->6Ky!;Jj8O9q4 z?A{%9+fd8cuivkYwk@`>VEoRIZwo|jzVOJ;N(b`pKr8DNU@5S~He^>WH*Tg8*!W5MR@pVZ;`@MYlM`UNMw4h&M={ovSY{;DHEJEi zWDc4rVa3H?3!kl%!wAnF7F^SlzGW8cG?KXAZCUZpaaxN2i}r9%Wtm1X#<+$)Z!c3L zW%CB|TXh>~@TnsNGd;>$i3!Z+5r-%)z@KXz_VAyngpgD z(+HmXO$ZUrWCh6yy@og?ka|X}ATkIeOwB*cBRdIbsmVH_)A?3D_t3p(g`<<=BLbuCD_7|q91SKced=h%W z30R|Y>K7I_!w{Ps_Tk+Fo-k3b$+i{0roy`W`h#`~-Kjp7)wgFzorjkYhwW1u?ztUA z4O;jWW0&a)ZxQO77no&+#WUMACY(QLkBiJ667I2r2C$8hiU$jIFRkEx3qq09x@%^R z?oXl=1pSv`JWDYZMdo4OBw1O-#v^yUI-}1$YcC$FH~4I~yfns3G->y8YiTdk9Y$(o z(by_UjB1bjRE_yLjZenr+6mgEg%Pu9I+kPnHrTf&&!Ptv9a7T=x+8? zePG)-QL?0hl~hqBm2s&akvahTMP2kQl@!k#El@>hc7xhd3D_p z^BI8oy`XfpPs4_ID`*fE{9nsCL%#NR8em=UqO{R;jX?TTA^4DSZ_g#Ez|2cdIu6~U zvn8{{xYz@IzE+VxLpFboJLez|KB-pG=Q+a;CWLV)G*=qJF5VjraQf1r;>v>hB&GUh zZT*Dj9W539T1@+^pnGB|(?b?d-r>}-cpUbZzNM{n7sWLj`f1;kn{Arwkvkkln|q&9 zZT>bMTPC$v_TCVr4i6G^yf2N?+>ULe z$eOi^&CO+L#E-@^Wt#DE8f{BsL)oKtx9mPCB9ujnv2Y|zncP+4mnCrGx8-0&6~VHQ zU-PNrOpuCRG*-d4?jESWb^K>#^+nP1kd%W9)ym+3T>7wkLh^1~2W~e1=%bdi=oenkqx5)dP8yLs zQfAf0r0T7;Wp#_M9WI{7ns6V-)aP(5q{1iv@?S(kOLW%g(buO#t`20o$h zcJFy+|7vQbC9!?Ki(q@j-`h(Uo$>&!fk#+Mq7ZItz;1cC}yOM7gJuku7UW@am~oSQk!#rvTDlk%nZGcud-~b9-g9g zf_^HDV{~A!d(j(aFFrL;{I^QDpe}w1woaEO^Ot3(`JA{53Wy4mJvKP zLo@z`cKo9;pG3(eaV&@Lr2BDAk#_8Dxb_4JE1tkO{tfy>ww-#GT8xs={nu6;OdXE{ zDv~WU*Wa^|=Ad+mh5!DS2KgjDIH4ijz0R0Zy1P!a2E|^O-3|x*p=b$CrLFs!)=)^_ zc?YT-abifrfwY#Q`KzPHRaX@OGVAwE3?kuNEX?7wN^M^w2 zsu9$`JXoL8W$OuVtGcU`=KUTWlRFVX<)`-GxUwVV&ummOtGtSUx&mjmWZ)u4h9@wo z<-TaS3e%T}_L%*{GRkadL(6`SHA^r(@%C7KW*CKB92EyCk*c%Y?__zb^wj1$ z`t;OA>U0meXi5%8=hc!26_Th@cESc}Seq7?b!?)3n=-1my52WW^$5{NDnb@r~zG^G#=2aR+T1Xsy%o#~c=pQ4ead2N4C(YgqCP{m&gmt%>gR4g4@zy?I2u?Bs4#iKyo8j7 zt2qx2i++;uDtSSdXw;=%iOeIEO{arAqHYU9%#(C*aUqecjdWO948GefRa~oB^GK){ zduuaK)s~eVa*HWv^l3s%jiftXu))<2|8$(aywOziR>E0*s!VUJ0eRTO-q9l`Pk0-G zQH7=4136{D0`W$*+di}OSTKP3Y+F7NA39+PT?41kx|o1;WDg;Y0QobvP_~5a+JfwJ zeG)wTAoLkrCkIFR3lto(JAFvQcW0in=0fiNVV*tpOMPjtW1Ada}inT7ow zvoc1>P3@V&L1-E5cG52QWk<}N2-DdDJ%U?;`Z!Cul#$FCmq6)U-M|bP@;qh(P4SLj zyjK)D!pHX^l*J0pO#~Vb@bM+JYVKjdu(kLl_GLzp^~U!GwlC1+2|okCby3(t;*x_Y zdiU)u?38ozKZ3m(M4zyp4|n4=&tCk|K?R>wtC``8Nih_Q5$9~2H~u!mq&$3+zI~9T zuFMh@JskkJpd!*PS2dxpi46cIOp#o_P zh#mrNglSXtA-6&?^>A0B12muYu_4+=W=d}9@SQa)1HhRHiP+t1=2)Ov+vuR3KH<%6 zzLyuAkHkS$S=#dv0(*WyTSlr7t&i`!a{f)uUP~W9j-fXX5{H94u1{0-&mU2^nv%82 zO`(k~u5);2f+bh;BS>D>+$bqUl1`o-G695@Jxwj^jC`T}Sp(8tzk5j$YZb zZK39x$qGTEX9t5yfwP3vcQkH^B6bz|;NpqdulE>#QThC0WIO_Ix92=;H(_H21|j>D zXM#**f4h`{=xmJXT9}vnylR|H8cGAPVYp{c;;68Ne+lg$dxT=LzQrKyakd~6;rhqn z*B}H*{Mu&P+Egv{VKyhie3$2c+lZ~&4ZB=9Dl%PWS*|&qElUgttL;`3P?snY=M50z zh7-=mR#|1bVynORC06RE)ZGDUi=Ej;lBi~EfsrF`7&l0n^j-?V2RcR=r;1g+#Y#^k z%NTC=VE-FhBr3LC+zaTktGWs{6A*h|m76sJ(Ni3bgF@@M1o?5%LR`Q$>|mdKHOTh< zul_w$=hOH$>*QS6DS?Wiz9!|-aw_O^$$R!X-)GP|xUyA^e)uR;?^BLThJ%Mj0fNJ` z_N_xM3Zveir$&D1M0;0E)DWYSbVI6wSnKqkR18sVdOjy5+MEHmk_84124UzA)H?33 z+>(T9+8d!AS};~CUGxAPs6t5vx36|T*kS5LnH&;W3`|xf@W`&3#$idRku3M2E5vAY zUwLm(lG?XAYL<)j*X(w?O={^hg@$#V*R^ltlI2M0n(x{I&KrU;`&8i6afuV8 zUd3{(8ab@p;@`P~3OGrxJyECrWW`53Y6k@PdwsZio^l{Kp@@=P7o)T?EYEeH8dmFV zv6`|TX>M(DYoK3iiR#!VxtS##8~nT!P|42v+db|V$4?~m8|dQYqla_hK-Qh|BwQ*b zl`OVW5g_AF{CoZUtE54cUp_$oww5I|)iQ#$DnQmk#X(bz%gIkA|F}&_N=4&Zaib}j zeO!F7f#G0Uxb3(MC*Nc;?+nBME$R{G8Op|LB7jM~<-mkLup)Z_xL<=@mTMX@8}sS! z7Te(iG3LVXkDcUvmF1c-ZNVb}97sVc(v=7)w@o~Ysc6n*t23M_s@4d!u>@qQF+_9h zG)|8cVDAW-4;5ts`j2~N^XAwHo%|;b)YS6k&rJEpW`|e|B7daRbfik!uFR813NOP8 zrZcEpra|Smw?$v3pO8$CbcMwqVo8jOo9dP9QPS*O9El#ufN(4+RjW-A|iDU;~x@vsmCuW zIHH3cq|6nFkA(@5ZA{rAbyUdrGSb64aVR=y#C^J`OvO<7xhf};-w-+vo-Eh*C zi$yDyX(1(1pe^+78;JeiaQNqvPxWY)`_;zIY1-0tq_ z{YJNBe=0QUeLgBZI-W5?Jh0_EUs#GhM1h!iNIakX2mi~j0isa>`zxAYbnn00%q`BC zFj^2dlowjz*y+?IzH(xwTdYdL7T;BvdF13H`FAlsi>TF{cz4%I7|_bBbDdtIt_J9c zc)_BlJ|>GmfedSAD~cC6#57YHhtd>EmbQc|PoQ@X+0wa(=an?ZTf^U9 z8pE{sc-pS3w8K>H6|qVZr1v;<07YxevW(c66~izA&I&sfk0F@C%w5~_mFB?_#p%{j zw|3^6@~#X8{1gT2abG^tfg7HwaVBxafpiL|8}0Hwoe?>~QhUWDghmrZByL8%ry6U8 zw8_syAu<{yZq4!xs3`tUFOU$rDdi%g-SLa|8l@UdftP#f!!SOLc7e$xNV7Z6Q4hW$ z*;aqWb_GXV&^!10yoE&i!>Up*Q*skdW*(%6^(P`l3pV)u`xdyHCK!p42Fb4a~2V1J2F)s$8^h^W-(`{?b z-HPca6_VwE`B9S1TAUbMg1TyPQ~d;Ov;q&=+#RKL%Bj7Z47?^EKjzbI(X>hOOJUh} zMz}(l0+`_|ZxW^I*mrDu_2R_z(tK0ts1&<+w?A0oW~t@76?FwP^=tHr%vTrpzwUVz zapX9jS!|lVF9e~V^@!O<^XqG)}3&O-j(xUNTVosz4%LVCaKik^p=lAcy z#Ds9(3;oKB=PFLy6ngQ6o_HwhIW=pIx>?4&ld{DQFWJ@sOmm3X{4Em9 zNa==#OwanM5O%jV;wWN z3ihfzc#mH_riLi-Gn?@J2HM}=`vkG7_m+UMGO9-G2K_eig;<-Gv8x+Bg_dV$PIpPN z&ab%yKk{_jf4E;ITT_c)cNaQtLJ?s405|0&Ai?acn<>PVv9YGBE!jdlFmmRLPuSG6 zLEVe_s84;;3?M}aU`idW@D2Ehk}Op3rK>tU8;L6`>*j=87dRI@Efif zuR37cQu9+(1UT*ib+KozAUqDZpmLp6T={k0mqf8t;&__tbq};HFJc z9B}Q&rkn=l%#??ps@+LAjMqmYMZ3#psM7B;MJLitPUZa9qXLT%tH=oa`|EE0ULhTi zw$9>C5|MA_gRQ!F+mLt>w|nIgdOsrLyI;jwHv?_TG8he(GJ)w+KYI)+ZNZe=KyZx`k6=kJqd0y5MqjcX6C0b>tsOP?U5Z<~K#? zDKAYC<9W=Y3Z|J}C3OJU8aIMm4 zZ6(IHE04eJb+w`6i@sB7xWD^itZkzX;XLQzVB_Vgjp^gxk=R{7C^@tY^(R};j}W~q zoVA!v_GUcRzo9JntEm^5oYhv*?~g1a)2^zlXWa~`dzkQbHwzA~`A_@N89Eui+roK@ z(A^B72=5OmUq5Aqy6s9S4zYu5Y_S798nL%pO0fW`5MNnCr&Y9>Gkz2r`8wG}nI79n zDAj>YLJTf=sGE0m)2sN$-HCvg@U91RG?09&Wq~%4I`?_NBxQ|JMPEmM!Bro&EwEccPUM^MNN;jvcJ1ziUk&TUn(zA;FyP zBqvMt(2W>Hlm6CI?`1+4^ivsBJ^7;B5zUIt>j;XFrr@TZ2vpn{8-RaIR)UQ1ZWQWlBIwOs_ z`hI!|K99Sq9%Ps=vaF(|wRU*qolHaUAE)rhF{z>-D&4~0t*m=l;>}##UR9ol65fO! zr{E#Kks5=US@;)tHQUMhLd?y4v_{^p+^`;yW?{oz7vAG#J&+}$y1R*aAbLTg1Ib?AufPgEXC$zeg#@nW^uh(B)(+WldThY*luFz3!}E^030z z4K_w`&YAK0#?K86@z0m-WO>o6D%Qqi@P9ArI-}Fr7JE4x)2nQ84}il?Wkw|iZk0e+ z*<3^Zd4i8aOWRaniuCtx2dvdD%a#B=)@7wBXukRvlf3x%5Q@oCdicXG|9b!yn+?yc8+mGxD+hN zzY`H@^{f4Y^3ww>HFDk|%mjlcKPo@x{{$8x|LQ^iBqblBjDPi@|6xhV2PowKgJ|O8 l^8b11#|NnEKe~wjKiz*iosaU@|1{l?${H1+J7*H$e*r|v2lM~{ literal 0 HcmV?d00001 diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 9a5250e8b25b7..e7a692b97d9af 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -31,10 +31,13 @@ #include #include #include +#include #include #include #include +#include "../src/RooFormulaUtils.h" + #include #include #include @@ -42,6 +45,8 @@ #include #include +#include +#include #include #include "gtest_wrapper.h" @@ -87,6 +92,33 @@ void randomizeParameters(const RooArgSet ¶meters) } } +/// Set or unset ROOFIT_FORMULA_BACKEND for the lifetime of the object, +/// resetting the backend's read-once cache on both ends. +class ScopedFormulaBackendEnv { +public: + ScopedFormulaBackendEnv(const char *value) + { + if (const char *old = std::getenv("ROOFIT_FORMULA_BACKEND")) + _old = old; + if (value) + setenv("ROOFIT_FORMULA_BACKEND", value, /*overwrite=*/1); + else + unsetenv("ROOFIT_FORMULA_BACKEND"); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + ~ScopedFormulaBackendEnv() + { + if (_old.empty()) + unsetenv("ROOFIT_FORMULA_BACKEND"); + else + setenv("ROOFIT_FORMULA_BACKEND", _old.c_str(), 1); + RooFormulaInternal::resetFormulaBackendForTesting(); + } + +private: + std::string _old; +}; + } // namespace using CreateNLLFunc = @@ -218,6 +250,90 @@ TEST_P(FactoryTest, NLLFit) EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, tol)); } +// Gradient agreement between the two RooFormula evaluation backends: for +// models containing RooFormulaVar/RooGenericPdf, the Clad gradient from the +// new codegen path that inlines the formula expression (AST backend) must +// agree with the gradient from the old path that calls the +// cling-JIT-compiled TFormula function by name (tformula backend), and with +// numerical differentiation of the reference NLL. +TEST(RooFuncWrapperFormula, GradientBackendAgreement) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + struct BackendResult { + double nll = 0.0; + std::vector grad; + std::vector numGrad; + }; + + // Restore the shared RooRandom generator state afterwards, so that this + // test does not perturb the data generated by the other tests in this + // binary (their fit-result comparisons are seed-sensitive). + TRandom3 savedRngState{*static_cast(RooRandom::randomGenerator())}; + + auto run = [](std::vector const &factoryExprs, const char *backend) { + ScopedFormulaBackendEnv env{backend}; + + RooWorkspace ws; + for (std::string const &expr : factoryExprs) { + ws.factory(expr); + } + + RooAbsPdf &model = *ws.pdf("model"); + RooArgSet observables{*ws.var("x")}; + + // Fixed seed: both backends must see bitwise-identical data. + RooRandom::randomGenerator()->SetSeed(1337); + std::unique_ptr data{model.generate(observables, 50)}; + + std::unique_ptr nllRef{model.createNLL(*data, RooFit::EvalBackend::Cpu())}; + std::unique_ptr nllFunc{model.createNLL(*data, RooFit::EvalBackend::Codegen())}; + + BackendResult res; + res.nll = nllFunc->getVal(); + + RooArgSet params; + nllFunc->getParameters(&observables, params); + res.grad.assign(params.size(), 0.0); + nllFunc->gradient(res.grad.data()); + for (auto *param : params) { + res.numGrad.push_back(getNumDerivative(*nllRef, static_cast(*param), observables)); + } + return res; + }; + + const std::vector> models{ + // Gaussian with RooFormulaVar parameters mixing operators and functions + {"expr::mu_shifted('mu + 0.1*sin(mu) + shift^2', {mu[1.0, -10, 10], shift[0.5, -10, 10]})", + "expr::sigma_scaled('sigma*(1 + 0.2*abs(alpha))', {sigma[2.0, 0.1, 10], alpha[0.3, -5, 5]})", + "Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"}, + // RooGenericPdf, normalized by numeric integration + {"EXPR::model('exp(-0.5*((x - m)/s)^2) + 0.1*(1 + 0.05*x)', {x[0, -10, 10], m[1.0, -10, 10], s[2.0, 0.1, " + "10]})"}}; + + for (auto const &factoryExprs : models) { + BackendResult ast = run(factoryExprs, "ast"); + BackendResult tformula = run(factoryExprs, "tformula"); + + ASSERT_EQ(ast.grad.size(), tformula.grad.size()); + ASSERT_FALSE(ast.grad.empty()); + + // Same data, same math: the two codegen paths must agree to near + // machine precision, with a small margin for Clad emitting differently + // associated derivative code across the JIT'd-function call boundary. + EXPECT_NEAR(ast.nll, tformula.nll, 1e-12 * std::abs(tformula.nll)); + for (std::size_t i = 0; i < ast.grad.size(); ++i) { + EXPECT_NEAR(ast.grad[i], tformula.grad[i], 1e-9 * std::max(1.0, std::abs(tformula.grad[i]))) + << "gradient component " << i; + // ... and both must be consistent with numerical differentiation. + EXPECT_NEAR(ast.grad[i], ast.numGrad[i], 1e-4 * std::max(1.0, std::abs(ast.numGrad[i]))) + << "gradient component " << i; + } + } + + *static_cast(RooRandom::randomGenerator()) = savedRngState; +} + /// Initial minimization that was not based on any other tutorial/test. FactoryTestParams param1{"Gaussian", [](RooWorkspace &ws) { From 9760a192a11086a2fa74a989c77d89178d03c9ec Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 28 Aug 2026 07:45:16 +0000 Subject: [PATCH 4/5] [RF] Add regression test for cling OOM with many formulas (#21052) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The out-of-memory failure in issue #21052 comes from formulas whose bodies differ only in their numeric literals (TRExFitter-style expression NormFactors like `(1-(SF_b*0.035)-(SF_c*0.138))/0.518` with per-region constants): each distinct-literal formula misses TFormula's JIT'd-function cache and costs two cling JIT compilations (the formula plus its validation clone) at roughly 130 kB and several milliseconds each, i.e. about 8 GB and minutes at the reported scale of 30000 formulas. The issue's attached reproducer does not show this, because its formula bodies all normalize to the same string and the cache absorbs them, so the new test bakes distinct literals into every formula. With the JIT-free expression backend, constructing and evaluating 30000 such RooFormulaVars takes about a second and tens of MB. The test asserts that every formula stays on the JIT-free backend (an empty uniqueFuncName(), which on the TFormula backend is always the name of the JIT-compiled function) plus generous order-of-magnitude wall-time and RSS-growth bounds. The codegen variant covers the second half of the issue: codegen used to force one TFormula JIT compilation per formula just to obtain a function name for the generated code to call. Now the expressions are inlined, so building and evaluating the likelihood of a model with 1000 distinct-literal formulas invokes cling once, for the whole squashed function: O(1) compilations per model instead of O(N) per formula. 🤖 Done with the help of AI --- roofit/roofitcore/test/CMakeLists.txt | 2 +- .../test/testRooFormulaEvaluator.cxx | 176 ++++++++++++++++++ 2 files changed, 177 insertions(+), 1 deletion(-) diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 23b3bd663f7f4..3ac8af02124ad 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -65,7 +65,7 @@ ROOT_ADD_GTEST(testRooDataSet testRooDataSet.cxx LIBRARIES Tree RooFitCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/dataSet_with_errors_6_26_10.root) ROOT_ADD_GTEST(testRooFormula testRooFormula.cxx LIBRARIES RooFitCore ROOT::TestSupport) ROOT_ADD_GTEST(testRooFormulaEvaluator testRooFormulaEvaluator.cxx testRooFormulaEvaluatorIO.cxx - LIBRARIES RooFitCore Hist MathCore + LIBRARIES RooFitCore RooFit Hist MathCore COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testRooFormulaEvaluator_legacy_ws.root) ROOT_ADD_GTEST(testProxiesAndCategories testProxiesAndCategories.cxx LIBRARIES RooFitCore diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx index 95b26f0f81cf3..e9f3591abac23 100644 --- a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -7,8 +7,12 @@ #include "../src/RooFormulaParser.h" #include "../src/RooExprEvaluator.h" +#include +#include #include +#include #include +#include #include #include @@ -19,15 +23,22 @@ #include +#include #include +#include #include #include #include #include +#include #include #include #include +#ifdef __linux__ +#include +#endif + namespace { /// Bitwise comparison; any two NaNs count as equal. @@ -1529,3 +1540,168 @@ TEST(RooFormulaEvaluator, EmittedCppRandomExpressions) } } } + +namespace { + +/// Resident set size of this process in kB, or -1 where unsupported. +long vmRSSkB() +{ +#ifdef __linux__ + std::ifstream in("/proc/self/status"); + std::string line; + while (std::getline(in, line)) { + if (line.rfind("VmRSS:", 0) == 0) + return std::atol(line.c_str() + 6); + } +#endif + return -1; +} + +/// A formula string with numeric literals that are distinct for each `i`, +/// mirroring the TRExFitter-style expression NormFactors reported in +/// https://github.com/root-project/root/issues/21052, e.g. +/// `(1-(SFb_pcbt90_20_40*0.035)-(SFc_pcbt85_20_40*0.138))/0.518`. +std::string distinctLiteralExpr(int i) +{ + char buf[64]; + const double k1 = 0.001 + 1e-6 * i; + const double k2 = 0.100 + 1e-6 * i; + std::snprintf(buf, sizeof(buf), "(1-(a*%.6f)-(b*%.6f))/0.518", k1, k2); + return buf; +} + +} // namespace + +// Regression test for https://github.com/root-project/root/issues/21052: +// constructing many RooFormulaVars whose formulas differ only in their +// numeric literals must not JIT-compile anything. On the old TFormula +// backend, each distinct-literal formula misses the JIT'd-function cache and +// costs two cling JIT compilations (the formula and its validation clone) at +// ~130 kB and ~6 ms each -- about 8 GB and several minutes for this N, which +// is the reported out-of-memory failure. Note that literals distinct per +// formula are essential: with identical formula bodies the JIT'd-function +// cache absorbs the repetition and the old path passes this test too. +TEST(RooFormulaEvaluator, ManyDistinctLiteralFormulas) +{ + // The default backend must handle this workload; shield the test from an + // ambient ROOFIT_FORMULA_BACKEND override. + ScopedBackendEnv env{nullptr}; + + constexpr int n = 30000; + + RooRealVar a("a", "a", 1.0, 0.1, 3.0); + RooRealVar b("b", "b", 1.0, 0.1, 3.0); + + const long rss0 = vmRSSkB(); + const auto t0 = std::chrono::steady_clock::now(); + + std::vector> fvs; + fvs.reserve(n); + for (int i = 0; i < n; ++i) { + const std::string name = "f_" + std::to_string(i); + fvs.emplace_back(std::make_unique(name.c_str(), distinctLiteralExpr(i).c_str(), RooArgList(a, b))); + RooFormulaVar &fv = *fvs.back(); + fv.getVal(); // force evaluation like a real fit setup would + // The whole batch must run on the JIT-free expression backend: a + // formula on the TFormula backend always reports the (non-empty) name + // of its cling-JIT-compiled function, so a non-empty name here means a + // JIT compilation happened for this formula. + ASSERT_TRUE(fv.getUniqueFuncName().empty()) << name; + } + + const auto t1 = std::chrono::steady_clock::now(); + const long rss1 = vmRSSkB(); + + // Spot-check values and that the expressions are emittable as inline C++ + // (the positive counterpart of the empty unique function name: exactly the + // JIT-free expression backend can emit C++). + auto varName = [](unsigned int i) { return "v[" + std::to_string(i) + "]"; }; + for (int i : {0, 12345, n - 1}) { + const double k1 = 0.001 + 1e-6 * i; + const double k2 = 0.100 + 1e-6 * i; + EXPECT_DOUBLE_EQ(fvs[i]->getVal(), (1. - k1 - k2) / 0.518) << i; + EXPECT_FALSE(fvs[i]->emitFormulaCpp(varName).empty()) << i; + } + + // Order-of-magnitude resource bounds, far above what the JIT-free path + // needs (measured: 1.1 s wall, 82 MB RSS growth) and far below what the + // per-formula JIT path would cost (minutes, gigabytes). + const double wallSeconds = std::chrono::duration(t1 - t0).count(); + EXPECT_LT(wallSeconds, 60.); + if (rss0 >= 0) { + EXPECT_LT((rss1 - rss0) / 1024., 400.) << "RSS growth in MB"; + } +} + +// The codegen counterpart of ManyDistinctLiteralFormulas, also for +// https://github.com/root-project/root/issues/21052: before the JIT-free +// backend, codegen forced one TFormula JIT compilation per RooFormulaVar just +// to obtain a function name for the generated code to call. Now the +// expressions are inlined into the generated code, so building and evaluating +// the likelihood of a model with many distinct-literal formulas involves no +// per-formula JIT. Cling is still invoked once, for the whole squashed +// likelihood function: the claim under test is O(1) compilations per model +// instead of O(N) per formula. +TEST(RooFormulaEvaluator, ManyDistinctLiteralFormulasCodegen) +{ + ScopedBackendEnv env{nullptr}; + + constexpr int n = 1000; + + RooRealVar x("x", "x", 0.0, -10.0, 10.0); + RooRealVar a("a", "a", 1.0, 0.1, 3.0); + RooRealVar b("b", "b", 1.0, 0.1, 3.0); + RooRealVar sigma("sigma", "sigma", 2.0, 0.1, 10.0); + + const long rss0 = vmRSSkB(); + const auto t0 = std::chrono::steady_clock::now(); + + // A Gaussian whose mean is the sum of n distinct-literal formula terms, + // each scaled such that the sum stays of order one. + std::vector> fvs; + RooArgList terms; + fvs.reserve(n); + for (int i = 0; i < n; ++i) { + const std::string name = "f_" + std::to_string(i); + const std::string expr = "(" + distinctLiteralExpr(i) + " - 0.9)/" + std::to_string(n) + ".0"; + fvs.emplace_back(std::make_unique(name.c_str(), expr.c_str(), RooArgList(a, b))); + terms.add(*fvs.back()); + } + RooAddition mean("mean", "mean", terms); + RooGaussian gauss("gauss", "gauss", x, mean, sigma); + + RooDataSet data("data", "data", x); + for (int i = 0; i < 20; ++i) { + x.setVal(-3.0 + 0.3 * i); + data.add(x); + } + + // CodegenNoGrad is the codegen backend without the (clad-dependent) + // gradient generation, which is not needed to test the compilation count. + std::unique_ptr nll{gauss.createNLL(data, RooFit::EvalBackend::CodegenNoGrad())}; + const double nllVal = nll->getVal(); + + const auto t1 = std::chrono::steady_clock::now(); + const long rss1 = vmRSSkB(); + + // Codegen must not have created any TFormula behind the scenes (it used + // to call getVal() and getUniqueFuncName() on each formula, forcing one + // JIT compilation per formula): every formula still has no JIT'd function + // name, i.e. it stayed on the JIT-free expression backend throughout. + for (auto const &fv : fvs) { + ASSERT_TRUE(fv->getUniqueFuncName().empty()) << fv->GetName(); + } + + // The generated code must agree with the regular CPU evaluation backend. + std::unique_ptr nllRef{gauss.createNLL(data, RooFit::EvalBackend::Cpu())}; + EXPECT_NEAR(nllVal, nllRef->getVal(), 1e-10 * std::abs(nllRef->getVal())); + + // Order-of-magnitude resource bounds (measured: 0.6 s wall, 40 MB RSS + // growth, dominated by the one cling compilation of the squashed + // likelihood function). + const double wallSeconds = std::chrono::duration(t1 - t0).count(); + EXPECT_LT(wallSeconds, 60.); + if (rss0 >= 0) { + EXPECT_LT((rss1 - rss0) / 1024., 400.) << "RSS growth in MB"; + } +} From 2107291142738efd0122be2a17378bc549b67fbb Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 28 Aug 2026 08:17:31 +0000 Subject: [PATCH 5/5] [RF] Vectorize batch evaluation of RooFormula AST programs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the scalar per-event loop in RooFormula::doEval() with chunked, vectorized evaluation of the JIT-free expression programs, closing the performance gap of the AST backend to the old cling-JIT path on long formulas and widening its lead on short ones. The instruction set of the postfix expression programs moves into RooBatchCompute (RooExprProgram.h); the same instruction vector still drives the scalar RooExprEvaluator::eval(). A new interface entry point, RooBatchCompute::computeExprProgram(), interprets a program with the loops interchanged: one instruction is applied across a chunk of RooBatchCompute::bufferSize (64) events at a time, with the value stack held as stack-allocated 64-double chunk buffers, so all intermediates stay L1-resident (the same chunking as the existing pdf kernels). The per-instruction loops are trivial elementwise operations that auto-vectorize in each of the GENERIC/SSE4.1/AVX/AVX2/AVX512 libraries, the interpreter dispatch amortizes over 64 events, per-event indirect calls disappear, and the scalar-vs-vector input broadcast decision hoists out of the per-event loop. Call instructions now carry their resolved function pointer; the exp/log/sin/cos/sqrt spelling families get dedicated opcodes so the batch interpreter can use the fast vectorizable VDT implementations (fast_exp etc.) where ROOT is built with VDT, exactly like the pdf kernels. RooFormula::doEval() now short-circuits the all-scalar case (every input span of size 1, the HistFactory expression-NormFactor shape of issue #21052) to a single scalar evaluation that is broadcast, routes expression programs through computeExprProgram(), and keeps the scalar per-event loop for the TFormula fallback backend and for programs deeper than the vector interpreter's fixed-size chunk stack. The evaluator remains free of mutable state, so concurrent doEval() calls stay safe. Numerical contract: without VDT every vectorized operation is the exact same double-precision operation per event as scalar evaluation, so batch results are bitwise identical (asserted in the new differential tests on such builds). With VDT, batch results can differ from scalar evaluation within RooBatchCompute's usual batch-vs-scalar tolerance (relative ~5e-14, see _toleranceCompareBatches in the vectorisedPDFs tests); the same class of difference already applies to every built-in pdf kernel. New differential tests compare the vectorized doEval() per event against scalar evaluation of the same program with mixed span sizes (vector observable, scalar parameters, an unused dependent with an empty span) across the real-world formula corpus, generated random expressions including ternary/comparison/logical operators, batch sizes around the 64-event chunking boundaries, a 10^6-event batch, NaN/Inf propagation, the all-scalar short-circuit, and the TFormula fallback. Benchmarks (AVX2 host, vdt=OFF build, RooFit::Evaluator driver, vs the scalar-AST doEval before this commit and the cling-JIT backend): a long ~30-instruction formula over 1M events drops from 166 ns/event to 64 ns/event, now faster than the JIT path (69.5); the short HistFactory shape x*p drops from 8.2 to 0.9 ns/event (JIT: 12.2), reaching 0.4-0.6 ns/event on cache-resident batches. An NLL over 1M events with a RooGenericPdf improves from 42 ms to 12.7 ms per evaluation (JIT backend: 23.5 ms). The vectorized path wins at every batch size in the sweep; identical output checksums against the JIT backend confirm the bitwise contract on VDT-less builds. 🤖 Done with the help of AI --- roofit/batchcompute/CMakeLists.txt | 1 + roofit/batchcompute/res/RooBatchCompute.h | 18 ++ roofit/batchcompute/res/RooExprProgram.h | 84 ++++++ roofit/batchcompute/src/Initialisation.cxx | 10 + roofit/batchcompute/src/RooBatchCompute.cxx | 192 ++++++++++++ roofit/roofitcore/src/RooExprEvaluator.cxx | 51 +++- roofit/roofitcore/src/RooExprEvaluator.h | 56 ++-- roofit/roofitcore/src/RooFormulaParser.cxx | 45 ++- roofit/roofitcore/src/RooFormulaUtils.cxx | 52 ++++ .../test/testRooFormulaEvaluator.cxx | 280 ++++++++++++++++++ 10 files changed, 743 insertions(+), 46 deletions(-) create mode 100644 roofit/batchcompute/res/RooExprProgram.h diff --git a/roofit/batchcompute/CMakeLists.txt b/roofit/batchcompute/CMakeLists.txt index b9b45820987ef..fc5f1e452a25d 100644 --- a/roofit/batchcompute/CMakeLists.txt +++ b/roofit/batchcompute/CMakeLists.txt @@ -27,6 +27,7 @@ if(NOT CMAKE_VERSION VERSION_LESS "3.23.0") BASE_DIRS res/ FILES res/RooBatchCompute.h + res/RooExprProgram.h ) target_sources( RooBatchCompute_GENERIC diff --git a/roofit/batchcompute/res/RooBatchCompute.h b/roofit/batchcompute/res/RooBatchCompute.h index bd7a7301b1d7a..bcc6e9247608d 100644 --- a/roofit/batchcompute/res/RooBatchCompute.h +++ b/roofit/batchcompute/res/RooBatchCompute.h @@ -13,6 +13,8 @@ #ifndef ROOFIT_BATCHCOMPUTE_ROOBATCHCOMPUTE_H #define ROOFIT_BATCHCOMPUTE_ROOBATCHCOMPUTE_H +#include "RooExprProgram.h" + #include #include //for R__EXTERN, needed for windows @@ -171,6 +173,15 @@ class RooBatchComputeInterface { virtual ~RooBatchComputeInterface() = default; virtual void compute(Config const &cfg, Computer, std::span output, VarSpan, ArgSpan) = 0; + /// Evaluate a postfix expression program (a formula compiled by RooFit's + /// JIT-free formula backend, see RooExprProgram.h) over a batch of events. + /// Input spans of size 1 are broadcast; `stackDepth` is the program's + /// maximum expression stack depth and must not exceed + /// maxExprProgramStackDepth. The default implementation throws; the CPU + /// backends implement it. + virtual void computeExprProgram(Config const &cfg, std::span code, unsigned int stackDepth, + std::span output, VarSpan vars); + virtual double reduceSum(Config const &cfg, InputArr input, size_t n) = 0; virtual ReduceNLLOutput reduceNLL(Config const &cfg, std::span probas, std::span weights, std::span offsetProbas) = 0; @@ -224,6 +235,13 @@ inline void compute(Config cfg, Computer comp, std::span output, compute(cfg, comp, output, VarSpan{vars.begin(), vars.end()}, extraArgs); } +inline void computeExprProgram(Config cfg, std::span code, unsigned int stackDepth, + std::span output, VarSpan vars) +{ + auto dispatch = cfg.useCuda() ? dispatchCUDA : dispatchCPU; + dispatch->computeExprProgram(cfg, code, stackDepth, output, vars); +} + inline double reduceSum(Config cfg, InputArr input, size_t n) { auto dispatch = cfg.useCuda() ? dispatchCUDA : dispatchCPU; diff --git a/roofit/batchcompute/res/RooExprProgram.h b/roofit/batchcompute/res/RooExprProgram.h new file mode 100644 index 0000000000000..6f700b478d5f6 --- /dev/null +++ b/roofit/batchcompute/res/RooExprProgram.h @@ -0,0 +1,84 @@ +/* + * Project: RooFit + * + * Copyright (c) 2026, CERN + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted according to the terms + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) + */ + +#ifndef ROOFIT_BATCHCOMPUTE_ROOEXPRPROGRAM_H +#define ROOFIT_BATCHCOMPUTE_ROOEXPRPROGRAM_H + +#include +#include + +namespace RooBatchCompute { + +/// Opcodes of the postfix expression programs compiled by RooFit's JIT-free +/// formula backend (see RooFormulaParser in RooFitCore). The same instruction +/// sequence drives both the scalar per-event evaluation in RooFitCore and the +/// chunked, vectorized batch evaluation in +/// RooBatchComputeInterface::computeExprProgram(). +enum class ExprOp : std::uint8_t { + Const, ///< push konst + Var, ///< push vars[arg] + Add, ///< a + b + Sub, ///< a - b + Mul, ///< a * b + Div, ///< a / b + Neg, ///< -a + Not, ///< !a (exactly 0.0 or 1.0) + LT, ///< a < b (exactly 0.0 or 1.0, likewise below) + LE, ///< a <= b + GT, ///< a > b + GE, ///< a >= b + EQ, ///< a == b + NE, ///< a != b + And, ///< a && b (no short-circuit: both operands are always evaluated) + Or, ///< a || b (no short-circuit) + Select, ///< c ? a : b (both branches are always evaluated) + Pow, ///< std::pow(a, b), from the `^`/`**` operator or pow() + Sq, ///< a * a, from TFormula's `expr^2` -> TMath::Sq(expr) rewrite + IntNorm, ///< a + 0.0: maps -0.0 to +0.0 where cling would have used integer arithmetic + // Unary calls whose semantics are exactly the corresponding std/libm + // function, split out from Call1 so that batch backends can substitute a + // fast vectorizable implementation (VDT, hardware sqrt). fn1 carries the + // exact scalar function, which is what per-event evaluation calls. + Exp, ///< std::exp(a) + Log, ///< std::log(a) + Sin, ///< std::sin(a) + Cos, ///< std::cos(a) + Sqrt, ///< std::sqrt(a) + Call1, ///< fn1(a) + Call2, ///< fn2(a, b) + Call3, ///< fn3(a, b, c) + Call4 ///< fn4(a, b, c, d) +}; + +/// One instruction of a postfix expression program. Call instructions carry +/// the resolved function pointer, so evaluation involves no lookup table; +/// `arg` additionally keeps the index into RooFitCore's function allow-list +/// (RooFormulaFunctions) that the call was resolved from, which C++ emission +/// uses to reproduce the exact spelling. +struct ExprInstr { + ExprOp op = ExprOp::Const; + std::uint32_t arg = 0; ///< Var: variable index; calls: function-table index + union { + double konst = 0.0; ///< Const + double (*fn1)(double); ///< Call1 and Exp...Sqrt + double (*fn2)(double, double); ///< Call2 + double (*fn3)(double, double, double); ///< Call3 + double (*fn4)(double, double, double, double); ///< Call4 + }; +}; + +/// Maximum expression stack depth accepted by computeExprProgram(), which +/// stack-allocates one bufferSize-sized chunk buffer per stack slot. Deeper +/// programs must be evaluated with the scalar per-event fallback. +constexpr std::uint32_t maxExprProgramStackDepth = 24; + +} // End namespace RooBatchCompute + +#endif diff --git a/roofit/batchcompute/src/Initialisation.cxx b/roofit/batchcompute/src/Initialisation.cxx index 391a81845b76b..df1dbf1f73710 100644 --- a/roofit/batchcompute/src/Initialisation.cxx +++ b/roofit/batchcompute/src/Initialisation.cxx @@ -46,6 +46,16 @@ bool &isInitialisedCpu() namespace RooBatchCompute { +/// Default implementation for backends that do not support evaluating +/// expression programs (currently the CUDA backend). RooFit only routes batch +/// formula evaluation through backends that do. +void RooBatchComputeInterface::computeExprProgram(Config const &, std::span, unsigned int, + std::span, VarSpan) +{ + throw std::runtime_error("computeExprProgram() is not implemented by the '" + architectureName() + + "' RooBatchCompute backend"); +} + /// Inspect hardware capabilities, and load the optimal library for RooFit computations. int initCPU() { diff --git a/roofit/batchcompute/src/RooBatchCompute.cxx b/roofit/batchcompute/src/RooBatchCompute.cxx index d98e0c828d6b5..ceb7d43f9e314 100644 --- a/roofit/batchcompute/src/RooBatchCompute.cxx +++ b/roofit/batchcompute/src/RooBatchCompute.cxx @@ -20,6 +20,7 @@ This file contains the code for cpu computations using the RooBatchCompute libra #include "RooBatchCompute.h" #include "RooNaNPacker.h" +#include "RooVDTHeaders.h" #include "Batches.h" #include @@ -31,6 +32,7 @@ This file contains the code for cpu computations using the RooBatchCompute libra #include #include +#include #include #include #include @@ -98,6 +100,8 @@ class RooBatchComputeClass : public RooBatchComputeInterface { }; void compute(Config const &, Computer computer, std::span output, VarSpan vars, ArgSpan extraArgs) override; + void computeExprProgram(Config const &, std::span code, unsigned int stackDepth, + std::span output, VarSpan vars) override; double reduceSum(Config const &, InputArr input, size_t n) override; ReduceNLLOutput reduceNLL(Config const &, std::span probas, std::span weights, std::span offsetProbas) override; @@ -234,6 +238,194 @@ void RooBatchComputeClass::compute(Config const &, Computer computer, std::span< namespace { +/// Apply a unary operation in place on a chunk at the top of the value stack. +template +inline void exprUnaryOp(double *__restrict a, std::size_t len, F f) +{ + for (std::size_t k = 0; k < len; ++k) { + a[k] = f(a[k]); + } +} + +/// Apply a binary operation on two chunks, storing the result in the first. +template +inline void exprBinaryOp(double *__restrict a, const double *__restrict b, std::size_t len, F f) +{ + for (std::size_t k = 0; k < len; ++k) { + a[k] = f(a[k], b[k]); + } +} + +} // namespace + +/** Evaluate a postfix expression program over a batch of events. + +The evaluation is chunked over bufferSize events, exactly like compute() and +the stack temporaries in ComputeFunctions.cxx, so that all intermediate value +buffers stay resident in L1 cache. Within a chunk, each instruction is applied +across the whole chunk: the per-instruction loops are trivial elementwise +operations that the compiler auto-vectorizes for the target architecture of +each RooBatchCompute library, and the interpreter dispatch cost is amortized +over bufferSize events. Scalar inputs (spans of size 1) are broadcast once per +chunk, hoisting the broadcast decision out of the per-event loop. + +Exp/Log/Sin/Cos use the fast vectorizable VDT implementations when ROOT is +built with VDT, exactly like the pdf compute kernels, in which case batch +results can differ from per-event scalar evaluation within the usual +RooBatchCompute batch-vs-scalar tolerance (relative ~5e-14). All other +operations apply the exact same double-precision operation per event that the +scalar evaluator applies, so without VDT the results are bitwise identical to +scalar evaluation. **/ +void RooBatchComputeClass::computeExprProgram(Config const &, std::span code, unsigned int stackDepth, + std::span output, VarSpan vars) +{ + if (stackDepth > maxExprProgramStackDepth) { + throw std::runtime_error("expression program exceeds the computeExprProgram() stack-depth limit"); + } + + double stack[maxExprProgramStackDepth][bufferSize]; + + const std::size_t nEvents = output.size(); + for (std::size_t begin = 0; begin < nEvents; begin += bufferSize) { + const std::size_t len = std::min(bufferSize, nEvents - begin); + std::size_t sp = 0; + for (ExprInstr const &ins : code) { + switch (ins.op) { + case ExprOp::Const: { + const double val = ins.konst; + double *__restrict out = stack[sp++]; + for (std::size_t k = 0; k < len; ++k) { + out[k] = val; + } + break; + } + case ExprOp::Var: { + std::span v = vars[ins.arg]; + double *__restrict out = stack[sp++]; + if (v.size() == 1) { + const double val = v[0]; + for (std::size_t k = 0; k < len; ++k) { + out[k] = val; + } + } else { + const double *__restrict in = v.data() + begin; + for (std::size_t k = 0; k < len; ++k) { + out[k] = in[k]; + } + } + break; + } + case ExprOp::Add: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a + b; }); + break; + case ExprOp::Sub: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a - b; }); + break; + case ExprOp::Mul: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a * b; }); + break; + case ExprOp::Div: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a / b; }); + break; + case ExprOp::Neg: exprUnaryOp(stack[sp - 1], len, [](double a) { return -a; }); break; + case ExprOp::Not: exprUnaryOp(stack[sp - 1], len, [](double a) { return a == 0.0 ? 1.0 : 0.0; }); break; + case ExprOp::LT: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a < b ? 1.0 : 0.0; }); + break; + case ExprOp::LE: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a <= b ? 1.0 : 0.0; }); + break; + case ExprOp::GT: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a > b ? 1.0 : 0.0; }); + break; + case ExprOp::GE: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a >= b ? 1.0 : 0.0; }); + break; + case ExprOp::EQ: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a == b ? 1.0 : 0.0; }); + break; + case ExprOp::NE: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return a != b ? 1.0 : 0.0; }); + break; + case ExprOp::And: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, + [](double a, double b) { return (a != 0.0 && b != 0.0) ? 1.0 : 0.0; }); + break; + case ExprOp::Or: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, + [](double a, double b) { return (a != 0.0 || b != 0.0) ? 1.0 : 0.0; }); + break; + case ExprOp::Select: { + sp -= 2; + double *__restrict c = stack[sp - 1]; + const double *__restrict a = stack[sp]; + const double *__restrict b = stack[sp + 1]; + for (std::size_t k = 0; k < len; ++k) { + c[k] = c[k] != 0.0 ? a[k] : b[k]; + } + break; + } + case ExprOp::Pow: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, [](double a, double b) { return std::pow(a, b); }); + break; + case ExprOp::Sq: exprUnaryOp(stack[sp - 1], len, [](double a) { return a * a; }); break; + case ExprOp::IntNorm: exprUnaryOp(stack[sp - 1], len, [](double a) { return a + 0.0; }); break; + case ExprOp::Exp: exprUnaryOp(stack[sp - 1], len, [](double a) { return fast_exp(a); }); break; + case ExprOp::Log: exprUnaryOp(stack[sp - 1], len, [](double a) { return fast_log(a); }); break; + case ExprOp::Sin: exprUnaryOp(stack[sp - 1], len, [](double a) { return fast_sin(a); }); break; + case ExprOp::Cos: exprUnaryOp(stack[sp - 1], len, [](double a) { return fast_cos(a); }); break; + case ExprOp::Sqrt: exprUnaryOp(stack[sp - 1], len, [](double a) { return std::sqrt(a); }); break; + case ExprOp::Call1: exprUnaryOp(stack[sp - 1], len, ins.fn1); break; + case ExprOp::Call2: + --sp; + exprBinaryOp(stack[sp - 1], stack[sp], len, ins.fn2); + break; + case ExprOp::Call3: { + sp -= 2; + double *__restrict a = stack[sp - 1]; + const double *__restrict b = stack[sp]; + const double *__restrict c = stack[sp + 1]; + for (std::size_t k = 0; k < len; ++k) { + a[k] = ins.fn3(a[k], b[k], c[k]); + } + break; + } + case ExprOp::Call4: { + sp -= 3; + double *__restrict a = stack[sp - 1]; + const double *__restrict b = stack[sp]; + const double *__restrict c = stack[sp + 1]; + const double *__restrict d = stack[sp + 2]; + for (std::size_t k = 0; k < len; ++k) { + a[k] = ins.fn4(a[k], b[k], c[k], d[k]); + } + break; + } + } + } + const double *__restrict res = stack[0]; + double *__restrict out = output.data() + begin; + for (std::size_t k = 0; k < len; ++k) { + out[k] = res[k]; + } + } +} + +namespace { + inline std::pair getLog(double prob, ReduceNLLOutput &out) { if (prob <= 0.0) { diff --git a/roofit/roofitcore/src/RooExprEvaluator.cxx b/roofit/roofitcore/src/RooExprEvaluator.cxx index 8044aca89f345..f3f4bf0797ceb 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.cxx +++ b/roofit/roofitcore/src/RooExprEvaluator.cxx @@ -92,6 +92,23 @@ Entry F4(const char *name, double (*fn)(double, double, double, double)) // (see TFormula::FillDefaults). All spellings below were cross-checked // against what TFormula accepts today; do not add a spelling TFormula // would reject (e.g. there is no std::sq or std::sign). +/// Tag the arity-1 entries with the given spellings with a dedicated opcode +/// (see Entry::op1). Every tagged spelling family calls exactly the same libm +/// function (the TMath:: versions are inline wrappers of the std ones), so the +/// scalar semantics carried by fn1 are unchanged; only the batch evaluation in +/// RooBatchCompute::computeExprProgram() treats the opcode specially. +void tagVectorizable(std::vector &table, std::initializer_list names, RooBatchCompute::ExprOp op) +{ + for (Entry &e : table) { + if (e.arity != 1) + continue; + for (const char *name : names) { + if (std::strcmp(e.name, name) == 0) + e.op1 = op; + } + } +} + std::vector makeTable() { auto castInt = +[](double x) { return static_cast(static_cast(x)); }; @@ -112,7 +129,7 @@ std::vector makeTable() auto sign = +[](double a, double b) { return TMath::Sign(a, b); }; auto signBit = +[](double x) { return std::signbit(x) ? 1.0 : 0.0; }; - return { + std::vector table{ // clang-format off // one-argument functions, libm/std spellings F1("sqrt", +[](double x) { return std::sqrt(x); }), @@ -222,6 +239,15 @@ std::vector makeTable() F4("TMath::Gaus", +[](double x, double m, double s, double n) { return TMath::Gaus(x, m, s, n != 0.0); }), // clang-format on }; + + using RooBatchCompute::ExprOp; + tagVectorizable(table, {"exp", "std::exp", "TMath::Exp"}, ExprOp::Exp); + tagVectorizable(table, {"log", "std::log", "TMath::Log"}, ExprOp::Log); + tagVectorizable(table, {"sin", "std::sin", "TMath::Sin"}, ExprOp::Sin); + tagVectorizable(table, {"cos", "std::cos", "TMath::Cos"}, ExprOp::Cos); + tagVectorizable(table, {"sqrt", "std::sqrt", "TMath::Sqrt"}, ExprOp::Sqrt); + + return table; } std::vector const &theTable() @@ -270,7 +296,6 @@ double RooExprEvaluator::eval(const double *vars) const double stack[kMaxStackDepth]; std::size_t sp = 0; - auto const *funcs = RooFormulaFunctions::table(); for (Instr const &ins : _program->code) { switch (ins.op) { @@ -336,18 +361,25 @@ double RooExprEvaluator::eval(const double *vars) const break; case Op::Sq: stack[sp - 1] *= stack[sp - 1]; break; case Op::IntNorm: stack[sp - 1] += 0.0; break; - case Op::Call1: stack[sp - 1] = funcs[ins.arg].fn1(stack[sp - 1]); break; + // The dedicated vectorizable opcodes carry the exact scalar function in + // fn1 (see Entry::op1), so scalar evaluation treats them like Call1. + case Op::Exp: + case Op::Log: + case Op::Sin: + case Op::Cos: + case Op::Sqrt: + case Op::Call1: stack[sp - 1] = ins.fn1(stack[sp - 1]); break; case Op::Call2: --sp; - stack[sp - 1] = funcs[ins.arg].fn2(stack[sp - 1], stack[sp]); + stack[sp - 1] = ins.fn2(stack[sp - 1], stack[sp]); break; case Op::Call3: sp -= 2; - stack[sp - 1] = funcs[ins.arg].fn3(stack[sp - 1], stack[sp], stack[sp + 1]); + stack[sp - 1] = ins.fn3(stack[sp - 1], stack[sp], stack[sp + 1]); break; case Op::Call4: sp -= 3; - stack[sp - 1] = funcs[ins.arg].fn4(stack[sp - 1], stack[sp], stack[sp + 1], stack[sp + 2]); + stack[sp - 1] = ins.fn4(stack[sp - 1], stack[sp], stack[sp + 1], stack[sp + 2]); break; } } @@ -458,6 +490,13 @@ std::string RooExprEvaluator::emitCpp(std::function c } case Op::Sq: stack.back() = "TMath::Sq(" + stack.back() + ")"; break; case Op::IntNorm: stack.back() = "(" + stack.back() + " + 0.0)"; break; + // The dedicated vectorizable opcodes keep the function-table index in + // `arg`, so they emit exactly like Call1 (with the original spelling). + case Op::Exp: + case Op::Log: + case Op::Sin: + case Op::Cos: + case Op::Sqrt: case Op::Call1: stack.back() = emissionName(funcs[ins.arg]) + "(" + stack.back() + ")"; break; case Op::Call2: { const std::string b = pop(); diff --git a/roofit/roofitcore/src/RooExprEvaluator.h b/roofit/roofitcore/src/RooExprEvaluator.h index cd1f1d9e424b4..5abb8a97e4b6c 100644 --- a/roofit/roofitcore/src/RooExprEvaluator.h +++ b/roofit/roofitcore/src/RooExprEvaluator.h @@ -15,6 +15,10 @@ #include "RooFormulaEvaluator.h" +#include "RooExprProgram.h" + +#include + #include #include #include @@ -53,6 +57,13 @@ struct Entry { const char *cppName = nullptr; std::uint8_t arity = 0; TypeRule rule = TypeRule::Double; + /// Opcode emitted for a call to this entry when arity == 1. Entries whose + /// semantics are exactly a std/libm function with a fast vectorizable + /// batch implementation (the exp/log/sin/cos/sqrt spelling families) get + /// the corresponding dedicated opcode instead of the generic Call1, so + /// that RooBatchCompute::computeExprProgram() can vectorize them. Scalar + /// evaluation and C++ emission treat those opcodes exactly like Call1. + RooBatchCompute::ExprOp op1 = RooBatchCompute::ExprOp::Call1; double (*fn0)() = nullptr; double (*fn1)(double) = nullptr; double (*fn2)(double, double) = nullptr; @@ -78,38 +89,12 @@ Entry const *find(std::string const &name, unsigned int nArgs, std::uint32_t &in /// call concurrently from multiple threads. class RooExprEvaluator final : public RooFormulaEvaluator { public: - enum class Op : std::uint8_t { - Const, ///< push konst - Var, ///< push vars[arg] - Add, ///< a + b - Sub, ///< a - b - Mul, ///< a * b - Div, ///< a / b - Neg, ///< -a - Not, ///< !a (exactly 0.0 or 1.0) - LT, ///< a < b (exactly 0.0 or 1.0, likewise below) - LE, ///< a <= b - GT, ///< a > b - GE, ///< a >= b - EQ, ///< a == b - NE, ///< a != b - And, ///< a && b (no short-circuit: both operands are always evaluated) - Or, ///< a || b (no short-circuit) - Select, ///< c ? a : b (both branches are always evaluated) - Pow, ///< std::pow(a, b), from the `^`/`**` operator or pow() - Sq, ///< a * a, from TFormula's `expr^2` -> TMath::Sq(expr) rewrite - IntNorm, ///< a + 0.0: maps -0.0 to +0.0 where cling would have used integer arithmetic - Call1, ///< RooFormulaFunctions::table()[arg].fn1 - Call2, ///< RooFormulaFunctions::table()[arg].fn2 - Call3, ///< RooFormulaFunctions::table()[arg].fn3 - Call4 ///< RooFormulaFunctions::table()[arg].fn4 - }; - - struct Instr { - Op op = Op::Const; - std::uint32_t arg = 0; - double konst = 0.0; - }; + /// The instruction set is shared with RooBatchCompute (see + /// RooExprProgram.h): the same instruction vector drives the scalar + /// per-event eval() here and the chunked, vectorized batch evaluation in + /// RooBatchCompute::computeExprProgram(). + using Op = RooBatchCompute::ExprOp; + using Instr = RooBatchCompute::ExprInstr; /// A compiled formula: immutable after construction and shared between all /// RooFormula instances with the same processed formula string. @@ -139,6 +124,13 @@ class RooExprEvaluator final : public RooFormulaEvaluator { std::string emitCpp(std::function const &varName) const override; + /// The compiled instruction sequence, for handing to + /// RooBatchCompute::computeExprProgram(). + std::span code() const { return {_program->code.data(), _program->code.size()}; } + + /// The program's maximum expression stack depth. + unsigned int stackDepth() const { return _program->stackDepth; } + private: std::shared_ptr _program; }; diff --git a/roofit/roofitcore/src/RooFormulaParser.cxx b/roofit/roofitcore/src/RooFormulaParser.cxx index 9ce2d3b650895..53a9a38cd9e64 100644 --- a/roofit/roofitcore/src/RooFormulaParser.cxx +++ b/roofit/roofitcore/src/RooFormulaParser.cxx @@ -464,7 +464,7 @@ class Parser { // over ExprOp with no default case, so that adding an opcode without // extending the accounting is a compiler warning (-Wswitch); the // static_assert additionally breaks the build when the enum grows. - static_assert(static_cast(Op::Call4) == 23, + static_assert(static_cast(Op::Call4) == 28, "ExprOp changed: update the stack-depth accounting switch in RooFormulaParser"); int depth = 0; int maxDepth = 0; @@ -476,6 +476,11 @@ class Parser { case Op::Not: case Op::Sq: case Op::IntNorm: + case Op::Exp: + case Op::Log: + case Op::Sin: + case Op::Cos: + case Op::Sqrt: case Op::Call1: break; case Op::Add: case Op::Sub: @@ -538,15 +543,39 @@ class Parser { Token const &peek() const { return _tokens[_pos]; } Token const &next() { return _tokens[_pos++]; } - void emit(Op op, std::uint32_t arg = 0, double konst = 0.0) + void emit(Op op, std::uint32_t arg = 0) { Instr ins; ins.op = op; ins.arg = arg; + _code.push_back(ins); + } + + void emitConst(double konst) + { + Instr ins; + ins.op = Op::Const; ins.konst = konst; _code.push_back(ins); } + /// Emit a call instruction. The resolved function pointer is stored in the + /// instruction itself (evaluation involves no table lookup); `index` keeps + /// the position in RooFormulaFunctions::table() for C++ emission. + void emitCall(Op op, std::uint32_t index, RooFormulaFunctions::Entry const &entry) + { + Instr ins; + ins.op = op; + ins.arg = index; + switch (op) { + case Op::Call2: ins.fn2 = entry.fn2; break; + case Op::Call3: ins.fn3 = entry.fn3; break; + case Op::Call4: ins.fn4 = entry.fn4; break; + default: ins.fn1 = entry.fn1; break; // Call1 and the vectorizable unary opcodes + } + _code.push_back(ins); + } + /// conditional-expression: right-associative, both branches always /// evaluated with the value of the active branch selected (Op::Select). bool parseTernary(ExprInfo &out) @@ -786,7 +815,7 @@ class Parser { switch (peek().kind) { case Tok::Number: { Token const &tok = next(); - emit(Op::Const, 0, tok.value); + emitConst(tok.value); out.type = tok.isInt ? ExprInfo::Type::Int : ExprInfo::Type::Double; out.isIntConst = tok.isInt; out.intConstValue = tok.intValue; @@ -883,12 +912,12 @@ class Parser { switch (nArgs) { case 0: // Zero-argument calls are the TMath constants: fold to a literal. - emit(Op::Const, 0, entry->fn0()); + emitConst(entry->fn0()); break; - case 1: emit(Op::Call1, index); break; - case 2: emit(Op::Call2, index); break; - case 3: emit(Op::Call3, index); break; - case 4: emit(Op::Call4, index); break; + case 1: emitCall(entry->op1, index, *entry); break; + case 2: emitCall(Op::Call2, index, *entry); break; + case 3: emitCall(Op::Call3, index, *entry); break; + case 4: emitCall(Op::Call4, index, *entry); break; } if (out.isIntegral()) emit(Op::IntNorm); // integer-valued calls cannot yield -0.0 in cling diff --git a/roofit/roofitcore/src/RooFormulaUtils.cxx b/roofit/roofitcore/src/RooFormulaUtils.cxx index 0e156e851f2ba..c50ccdb0a242f 100644 --- a/roofit/roofitcore/src/RooFormulaUtils.cxx +++ b/roofit/roofitcore/src/RooFormulaUtils.cxx @@ -33,6 +33,10 @@ the parser does not support silently fall back to the traditional TFormula environment variable `ROOFIT_FORMULA_BACKEND` overrides this: `tformula` always uses the TFormula backend, and `ast` disables the fallback, turning unsupported expressions into hard errors (useful for testing). + +Batch evaluation (doEvalFormula()) of expressions on the built-in backend is +chunk-vectorized via the RooBatchCompute library, see doEvalFormula() for the +numerical implications. **/ #include "RooFormulaUtils.h" @@ -44,10 +48,12 @@ unsupported expressions into hard errors (useful for testing). #include "RooCurve.h" #include "RooFitImplHelpers.h" #include "RooMsgService.h" +#include "RooBatchCompute.h" #include "RooExprEvaluator.h" #include "RooFormulaParser.h" #include "RooTFormulaEvaluator.h" +#include #include #include #include @@ -534,21 +540,67 @@ RooFormulaUtils::evalFormula(RooFormulaEvaluator const &evaluator, RooAbsCollect //////////////////////////////////////////////////////////////////////////////// /// Evaluate a formula for a batch of input values from the evaluation context, /// with `x[i]` taking the values of the i-th variable in `actualVars`. +/// +/// If every input is a single value (e.g. a formula of parameters only, the +/// HistFactory expression-NormFactor shape), the formula is evaluated once and +/// the result is broadcast. Otherwise, formulas on the JIT-free expression +/// backend are evaluated with RooBatchCompute::computeExprProgram(), which +/// applies one instruction across a chunk of RooBatchCompute::bufferSize +/// events at a time in vectorizable elementwise loops. When ROOT is built with +/// VDT, that path evaluates exp/log/sin/cos with the same fast vectorizable +/// implementations as the RooBatchCompute pdf kernels, so batch results can +/// differ from per-event scalar evaluation within the usual RooBatchCompute +/// batch-vs-scalar tolerance (relative ~5e-14, see the vectorisedPDFs tests); +/// without VDT the batch results are bitwise identical to scalar evaluation. +/// The TFormula fallback backend evaluates with a scalar per-event loop as +/// before. void RooFormulaUtils::doEvalFormula(RooFormulaEvaluator const &evaluator, RooArgList const &actualVars, RooFit::EvalContext &ctx) { std::span output = ctx.output(); + // Every x[i] has an input span, because the evaluation engine refers to + // the (pruned) list of actual dependents directly. const std::size_t nPars = actualVars.size(); // Note: emplace_back() instead of assignment into a pre-sized vector, // because the custom std::span backport for C++ < 20 in ROOT/span.hxx is // not move-assignable. std::vector> inputSpans; inputSpans.reserve(nPars); + bool allScalar = true; for (std::size_t i = 0; i < nPars; ++i) { inputSpans.emplace_back(ctx.at(static_cast(&actualVars[i]))); + allScalar &= inputSpans.back().size() <= 1; + } + + // All inputs are single values: evaluate once and broadcast. + if (allScalar) { + std::vector pars(nPars); + for (std::size_t j = 0; j < nPars; j++) { + if (!inputSpans[j].empty()) { + pars[j] = inputSpans[j][0]; + } + } + std::fill(output.begin(), output.end(), evaluator.eval(pars.data())); + return; + } + + // Chunked, vectorized evaluation of JIT-free expression programs. + if (auto *expr = dynamic_cast(&evaluator)) { + if (expr->stackDepth() <= RooBatchCompute::maxExprProgramStackDepth) { + // Load the RooBatchCompute CPU dispatch if no RooFit::Evaluator has + // done so yet: doEvalFormula() can also be called directly, and the + // dispatch pointer is null until the library is loaded (a cheap + // no-op once initialized). + RooBatchCompute::initCPU(); + RooBatchCompute::computeExprProgram({}, expr->code(), expr->stackDepth(), output, + {inputSpans.data(), inputSpans.size()}); + return; + } } + // Scalar per-event loop: the TFormula fallback backend, and expression + // programs too deep for the vector interpreter's fixed-size chunk stack. std::vector pars(nPars); for (std::size_t i = 0; i < output.size(); ++i) { for (std::size_t j = 0; j < nPars; ++j) { diff --git a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx index e9f3591abac23..eef541fd509b5 100644 --- a/roofit/roofitcore/test/testRooFormulaEvaluator.cxx +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -16,6 +17,7 @@ #include #include +#include // for R__HAS_VDT #include #include #include @@ -1543,6 +1545,284 @@ TEST(RooFormulaEvaluator, EmittedCppRandomExpressions) namespace { +// The vectorized doEval() path evaluates exp/log/sin/cos with the fast VDT +// implementations when ROOT is built with VDT (like all RooBatchCompute pdf +// kernels), in which case batch and per-event scalar results agree within +// RooBatchCompute's own batch-vs-scalar tolerance (_toleranceCompareBatches +// in roofit/test/vectorisedPDFs/VectorisedPDFTests.h). Without VDT, every +// vectorized operation is the exact same double-precision operation the +// scalar evaluator applies, so agreement must be bitwise. +bool batchAgrees(double batch, double ref) +{ + if (sameBits(batch, ref)) { + return true; + } +#ifdef R__HAS_VDT + if (std::isnan(batch) || std::isnan(ref) || std::isinf(batch) || std::isinf(ref)) { + return false; + } + return std::abs(batch - ref) <= 5e-14 * std::max(1.0, std::abs(ref)); +#else + return false; +#endif +} + +/// Batch-evaluate `processedExpr` (in the x[i] dialect) through +/// RooFit::Evaluator, exercising RooFormula::doEval() with mixed span sizes: +/// x[0] gets the vector input `xData` (span of size N), all other x[i] are +/// scalar parameters (spans of size 1) with value scalarVals[i], and one +/// trailing unused dependent exercises the empty-span handling. +std::vector batchEvalFormula(std::string const &processedExpr, std::vector const &xData, + std::vector const &scalarVals) +{ + const std::size_t nVars = scalarVals.size(); + RooArgList vars; + std::vector> owned; + for (std::size_t i = 0; i <= nVars; ++i) { // one extra, unused dependent + const std::string name = "v" + std::to_string(i); + const double val = i == 0 ? (xData.empty() ? 1.0 : xData[0]) : (i < nVars ? scalarVals[i] : 0.5); + owned.emplace_back(std::make_unique(name.c_str(), name.c_str(), val, -1e300, 1e300)); + vars.add(*owned.back()); + } + RooFormulaVar f("f", processedExpr.c_str(), vars); + RooFit::Evaluator ev(f); + ev.setInput("v0", {xData.data(), xData.size()}, false); + std::span out = ev.run(); + return {out.begin(), out.end()}; +} + +/// Per-event scalar reference for the same inputs, through the same compiled +/// program that the vectorized path executes. +std::vector scalarRefFormula(std::string const &processedExpr, std::vector const &xData, + std::vector const &scalarVals) +{ + auto prog = RooFormulaParser::compile(processedExpr, scalarVals.size() + 1); + if (!prog) { + ADD_FAILURE() << "expression unexpectedly failed to parse: " << processedExpr; + return {}; + } + RooExprEvaluator ev{prog}; + std::vector pars(scalarVals.size() + 1); + for (std::size_t i = 1; i < scalarVals.size(); ++i) { + pars[i] = scalarVals[i]; + } + std::vector out(xData.size()); + for (std::size_t i = 0; i < xData.size(); ++i) { + pars[0] = xData[i]; + out[i] = ev.eval(pars.data()); + } + return out; +} + +/// Compare a batch output against the per-event reference. An output of size +/// 1 means the evaluator collapsed the case to a single value (x[0] unused or +/// a size-1 input span): compare only the first reference value then. +void expectBatchMatches(std::vector const &out, std::vector const &ref, std::string const &what) +{ + ASSERT_FALSE(ref.empty()) << what; + if (out.size() == 1) { + EXPECT_TRUE(batchAgrees(out[0], ref[0])) + << what << "\n batch = " << std::hexfloat << out[0] << " scalar = " << ref[0] << std::defaultfloat; + return; + } + ASSERT_EQ(out.size(), ref.size()) << what; + for (std::size_t i = 0; i < out.size(); ++i) { + ASSERT_TRUE(batchAgrees(out[i], ref[i])) << what << "\n event " << i << ": batch = " << std::hexfloat << out[i] + << " scalar = " << ref[i] << std::defaultfloat; + } +} + +} // namespace + +// Differential test of the vectorized doEval() against per-event scalar +// evaluation on batch sizes around the bufferSize=64 chunking boundaries, +// with an expression covering ternary, comparisons, logical operators and the +// vectorizable functions, and with input values that produce NaN (log of a +// negative number, sqrt of a negative number) and Inf (division by zero), to +// check that special values propagate identically. +TEST(RooFormulaEvaluator, VectorizedDoEvalEdgeSizes) +{ + ScopedBackendEnv env{nullptr}; + + const std::string expr = "x[0]*x[1] + sin(x[0])*cos(x[2]) + (x[0] > 0.5 ? log(x[0] - 1.0) : -x[0])" + " + sqrt(x[0] - 2.0) + (x[0] != 0.0 && x[1] > 0.0) + 1.0/x[0] + exp(-x[0])"; + const std::vector scalarVals{0.0, 1.5, 0.7}; + + std::mt19937 rng{20260828u}; + std::uniform_real_distribution dist{-3.0, 8.0}; + + for (std::size_t n : {1u, 2u, 63u, 64u, 65u, 127u, 128u, 1000u}) { + std::vector xData(n); + for (double &v : xData) { + v = dist(rng); + } + // special values: division by zero -> Inf, negatives -> NaN from log/sqrt + if (n > 2) { + xData[0] = 0.0; + xData[1] = -1.0; + xData[2] = 1.75; // log(0.75) finite, sqrt(-0.25) NaN + } + auto out = batchEvalFormula(expr, xData, scalarVals); + auto ref = scalarRefFormula(expr, xData, scalarVals); + expectBatchMatches(out, ref, "n = " + std::to_string(n)); + } +} + +// Differential doEval() over the real-world corpus: every entry the JIT-free +// parser accepts is evaluated through the batch path (x[0] vectorized, other +// variables scalar, plus an unused dependent) and compared per event against +// scalar evaluation of the same program. +TEST(RooFormulaEvaluator, VectorizedDoEvalCorpus) +{ + ScopedBackendEnv env{nullptr}; + + constexpr std::size_t nEvents = 197; // 3 full chunks plus a remainder + + std::mt19937 rng{555u}; + std::uniform_real_distribution dist{0.1, 3.0}; + + for (const char *entry : kCorpus) { + int nVars = 0; + const std::string processed = normalizeCorpusEntry(entry, nVars); + auto prog = RooFormulaParser::compile(processed, std::max(nVars, 1)); + if (!prog) { + continue; // TFormula-fallback expressions are not vectorized + } + + std::vector xData(nEvents); + for (double &v : xData) { + v = dist(rng); + } + xData[0] = 0.0; + xData[1] = -1.5; + + std::vector scalarVals(std::max(nVars, 1)); + for (std::size_t i = 1; i < scalarVals.size(); ++i) { + scalarVals[i] = dist(rng); + } + + auto out = batchEvalFormula(processed, xData, scalarVals); + auto ref = scalarRefFormula(processed, xData, scalarVals); + expectBatchMatches(out, ref, std::string{"corpus entry: "} + entry + "\n processed: " + processed); + } +} + +// Differential doEval() on randomly generated expressions (same generator and +// seed as RandomExpressionDifferential), covering ternary/comparison/logical +// operators and the whole function allow-list in random combinations. +TEST(RooFormulaEvaluator, VectorizedDoEvalRandomExpressions) +{ + ScopedBackendEnv env{nullptr}; + + constexpr int nExprs = 150; + constexpr std::size_t nEvents = 130; // two full chunks plus a remainder + + RandomExprGenerator gen{kRandomExprSeed, kRandomExprVars}; + std::mt19937 inputRng{424242u}; + std::uniform_real_distribution dist{-5.0, 5.0}; + + for (int iExpr = 0; iExpr < nExprs; ++iExpr) { + const RandomExpr expr = gen.gen(4); + + std::vector xData(nEvents); + for (double &v : xData) { + v = dist(inputRng); + } + xData[0] = 0.0; + xData[1] = 1e300; + xData[2] = -1e-300; + + std::vector scalarVals(kRandomExprVars); + for (std::size_t i = 1; i < scalarVals.size(); ++i) { + scalarVals[i] = dist(inputRng); + } + + auto out = batchEvalFormula(expr.text, xData, scalarVals); + auto ref = scalarRefFormula(expr.text, xData, scalarVals); + expectBatchMatches(out, ref, "random expression: " + expr.text); + } +} + +// A large batch through the vectorized path: the ~30-instruction expression +// from the Phase 5 benchmarks over 10^6 events, compared per event. +TEST(RooFormulaEvaluator, VectorizedDoEvalLargeBatch) +{ + ScopedBackendEnv env{nullptr}; + + const std::string expr = "0.5*exp(-0.5*(x[0]-x[1])*(x[0]-x[1])/(x[2]*x[2])) + 0.3*sin(0.5*x[0]+x[1])*cos(x[0]*x[2])" + " + 0.2/(1.0+x[0]*x[0]) + sqrt(abs(x[0]*x[1])+1.0) + 0.1*log(1.0+exp(-x[0]))"; + const std::vector scalarVals{0.0, 1.0, 2.0}; + + constexpr std::size_t nEvents = 1000000; + std::vector xData(nEvents); + std::mt19937 rng{31415u}; + std::uniform_real_distribution dist{0.0, 10.0}; + for (double &v : xData) { + v = dist(rng); + } + + auto out = batchEvalFormula(expr, xData, scalarVals); + auto ref = scalarRefFormula(expr, xData, scalarVals); + expectBatchMatches(out, ref, "large batch"); +} + +// The all-scalar case (every input span has size 1, e.g. a formula of +// parameters only -- the HistFactory NormFactor shape of issue #21052) must +// short-circuit to a single scalar evaluation, bitwise identical to eval(), +// with or without VDT. +TEST(RooFormulaEvaluator, VectorizedDoEvalAllScalar) +{ + ScopedBackendEnv env{nullptr}; + + RooRealVar a("a", "a", 1.3, 0.0, 10.0); + RooRealVar b("b", "b", 0.7, 0.0, 10.0); + RooRealVar unused("unused", "unused", 2.0, 0.0, 10.0); + RooFormulaVar f("f", "exp(-a*b) + sin(a)/b", {a, b, unused}); + + const double refVal = f.getVal(); + RooFit::Evaluator ev(f); + std::span out = ev.run(); + ASSERT_EQ(out.size(), 1u); + EXPECT_TRUE(sameBits(out[0], refVal)) << std::hexfloat << out[0] << " vs " << refVal << std::defaultfloat; +} + +// The TFormula fallback backend keeps its scalar per-event loop in doEval(); +// its batch results must stay bitwise identical to scalar evaluation. +TEST(RooFormulaEvaluator, TFormulaBackendDoEval) +{ + ScopedBackendEnv env{"tformula"}; + + RooRealVar x("x", "x", 5.0, 0.0, 10.0); + RooRealVar p("p", "p", 1.5, 0.1, 3.0); + RooRealVar q("q", "q", 0.7, 0.1, 3.0); // unused by the formula + RooFormulaVar f("f", "x*p + sin(x) + (x > 5.0 ? log(x) : -x)", {x, p, q}); + + std::vector xData(100); + std::mt19937 rng{777u}; + std::uniform_real_distribution dist{0.0, 10.0}; + for (double &v : xData) { + v = dist(rng); + } + + RooFit::Evaluator ev(f); + ev.setInput("x", {xData.data(), xData.size()}, false); + std::span out = ev.run(); + // The TFormula backend was used: there is a JIT-compiled function. + EXPECT_FALSE(f.getUniqueFuncName().empty()); + ASSERT_EQ(out.size(), xData.size()); + // Reference through the JIT-free scalar program, which evaluates bitwise + // identically to the JIT-compiled TFormula (the Phase 3 contract). A + // reference expression inlined here would not be reliable: this test file + // is compiled with FMA contraction enabled. + auto ref = scalarRefFormula("x[0]*x[1] + sin(x[0]) + (x[0] > 5.0 ? log(x[0]) : -x[0])", xData, {0.0, 1.5}); + ASSERT_EQ(ref.size(), xData.size()); + for (std::size_t i = 0; i < xData.size(); ++i) { + EXPECT_TRUE(sameBits(out[i], ref[i])) << "event " << i; + } +} + +namespace { + /// Resident set size of this process in kB, or -1 where unsupported. long vmRSSkB() {