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/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/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..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 @@ -83,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; diff --git a/roofit/roofitcore/inc/RooGenericPdf.h b/roofit/roofitcore/inc/RooGenericPdf.h index b76158fb85979..9a6598a4dd924 100644 --- a/roofit/roofitcore/inc/RooGenericPdf.h +++ b/roofit/roofitcore/inc/RooGenericPdf.h @@ -20,6 +20,7 @@ #include "RooListProxy.h" #include "RooAbsBinning.h" +#include #include #include #include @@ -65,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; diff --git a/roofit/roofitcore/src/RooExprEvaluator.cxx b/roofit/roofitcore/src/RooExprEvaluator.cxx new file mode 100644 index 0000000000000..f3f4bf0797ceb --- /dev/null +++ b/roofit/roofitcore/src/RooExprEvaluator.cxx @@ -0,0 +1,525 @@ +/// \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 +#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, 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, + const char *cppName = nullptr) +{ + Entry e; + e.name = name; + e.cppName = cppName; + 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). +/// 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)); }; + 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 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; }; + + std::vector table{ + // 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); }), + // 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); }), + 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::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); }), + 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), + // TFormula shortcut for TMath::Sign + 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(); }), + 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 + }; + + 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() +{ + 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; + + 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; + // 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] = ins.fn2(stack[sp - 1], stack[sp]); + break; + case Op::Call3: + sp -= 2; + stack[sp - 1] = ins.fn3(stack[sp - 1], stack[sp], stack[sp + 1]); + break; + case Op::Call4: + sp -= 3; + stack[sp - 1] = ins.fn4(stack[sp - 1], stack[sp], stack[sp + 1], stack[sp + 2]); + break; + } + } + + 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; + // 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(); + 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 new file mode 100644 index 0000000000000..5abb8a97e4b6c --- /dev/null +++ b/roofit/roofitcore/src/RooExprEvaluator.h @@ -0,0 +1,140 @@ +/// \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 "RooExprProgram.h" + +#include + +#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 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, ///< 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 { + 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; + /// 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; + 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: + /// 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. + 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; } + + bool canEmitCpp() const override { return true; } + + 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; +}; + +#endif + +/// \endcond 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/RooFormulaParser.cxx b/roofit/roofitcore/src/RooFormulaParser.cxx new file mode 100644 index 0000000000000..53a9a38cd9e64 --- /dev/null +++ b/roofit/roofitcore/src/RooFormulaParser.cxx @@ -0,0 +1,987 @@ +/// \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 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 + * floating-point exceptions are trapped. + */ + +#include "RooFormulaParser.h" + +#include +#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 + 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 +}; + +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"); + // 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) + 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 '+': + // 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; + 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++ 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 { + 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) + { + _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) == 28, + "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::Exp: + case Op::Log: + case Op::Sin: + case Op::Cos: + case Op::Sqrt: + 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; + }; + + 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()) + _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) + { + 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) + { + 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 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; + } + + /// Precedence climbing over gBinaryOps. + bool parseBinary(int minPrec, ExprInfo &out) + { + 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)) + return false; + switch (info->op) { + case Op::Add: + case Op::Sub: + case Op::Mul: + // 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.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.isIntegral() && rhs.isIntegral()) + return fail("integer division has truncating semantics in TFormula/cling"); + out.type = ExprInfo::Type::Double; + out.isIntConst = false; + emit(info->op); + break; + default: + // comparisons and logical operators: C++ result type is bool + out.type = ExprInfo::Type::Bool; + out.isIntConst = false; + 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(); + 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.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.type = ExprInfo::Type::Bool; + out.isIntConst = false; + 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.type = ExprInfo::Type::Double; // pow() and TMath::Sq(Double_t) return double + out.isIntConst = false; + 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; + } + // 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; + // 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.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; + } + + /// 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(); + emitConst(tok.value); + out.type = tok.isInt ? ExprInfo::Type::Int : ExprInfo::Type::Double; + out.isIntConst = tok.isInt; + out.intConstValue = tok.intValue; + 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.type = ExprInfo::Type::Double; + out.isIntConst = 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; + using Type = ExprInfo::Type; + switch (entry->rule) { + 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 (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: + // Zero-argument calls are the TMath constants: fold to a literal. + emitConst(entry->fn0()); + 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 + 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..c50ccdb0a242f 100644 --- a/roofit/roofitcore/src/RooFormulaUtils.cxx +++ b/roofit/roofitcore/src/RooFormulaUtils.cxx @@ -24,6 +24,19 @@ 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). + +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" @@ -35,14 +48,18 @@ RooFit::DEBUG message level for the RooFit::InputArguments topic. #include "RooCurve.h" #include "RooFitImplHelpers.h" #include "RooMsgService.h" +#include "RooBatchCompute.h" +#include "RooExprEvaluator.h" +#include "RooFormulaParser.h" #include "RooTFormulaEvaluator.h" -#include "TFormula.h" - +#include #include #include +#include #include #include +#include #include #include @@ -50,6 +67,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 +404,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 +423,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); } @@ -385,15 +505,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; } @@ -422,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/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..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" @@ -326,9 +324,40 @@ double RooFormulaVar::defaultErrorLevel() const return 1.0 ; } +//////////////////////////////////////////////////////////////////////////////// +/// 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 { - return evaluator().getTFormula()->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 549c293995f1b..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; @@ -227,7 +225,38 @@ void RooGenericPdf::writeToStream(ostream& os, bool compact) const } } +//////////////////////////////////////////////////////////////////////////////// +/// 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 { - return evaluator().getTFormula()->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/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index ac9d5bb73b7f0..3ac8af02124ad 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -64,6 +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 testRooFormulaEvaluatorIO.cxx + LIBRARIES RooFitCore RooFit 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 new file mode 100644 index 0000000000000..eef541fd509b5 --- /dev/null +++ b/roofit/roofitcore/test/testRooFormulaEvaluator.cxx @@ -0,0 +1,1987 @@ +// 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 // for R__HAS_VDT +#include +#include +#include +#include + +#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. +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)")); + // `^` 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")); + 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)); +} + +// 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. +// 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))); + + // 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) +{ + 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))); + // 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"}; + 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); + } +} + +// 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 +// 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"; + } + 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); +} + +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; +} + +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; + } + } +} + +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() +{ +#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"; + } +} 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 0000000000000..9b17c3d055d70 Binary files /dev/null and b/roofit/roofitcore/test/testRooFormulaEvaluator_legacy_ws.root differ 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) {