[RF] Implement RooFormulaVar and RooGenericPdf evalutation path without TFormula

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

roofit/batchcompute/res/RooBatchCompute.h

@@ -13,6 +13,8 @@
 #ifndef ROOFIT_BATCHCOMPUTE_ROOBATCHCOMPUTE_H
 #define ROOFIT_BATCHCOMPUTE_ROOBATCHCOMPUTE_H
 
+#include "RooExprProgram.h"
+
 #include <ROOT/RSpan.hxx>
 
 #include <DllImport.h> //for R__EXTERN, needed for windows
@@ -171,6 +173,15 @@ class RooBatchComputeInterface {
    virtual ~RooBatchComputeInterface() = default;
    virtual void compute(Config const &cfg, Computer, std::span<double> 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<const ExprInstr> code, unsigned int stackDepth,
+                                   std::span<double> output, VarSpan vars);
+
    virtual double reduceSum(Config const &cfg, InputArr input, size_t n) = 0;
    virtual ReduceNLLOutput reduceNLL(Config const &cfg, std::span<const double> probas, std::span<const double> weights,
                                      std::span<const double> offsetProbas) = 0;
@@ -224,6 +235,13 @@ inline void compute(Config cfg, Computer comp, std::span<double> output,
    compute(cfg, comp, output, VarSpan{vars.begin(), vars.end()}, extraArgs);
 }
 
+inline void computeExprProgram(Config cfg, std::span<const ExprInstr> code, unsigned int stackDepth,
+                               std::span<double> 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;

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 <cstddef>
+#include <cstdint>
+
+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

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<const ExprInstr>, unsigned int,
+                                                  std::span<double>, 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()
 {

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 <ROOT/RConfig.hxx>
@@ -31,6 +32,7 @@ This file contains the code for cpu computations using the RooBatchCompute libra
 #include <Math/Util.h>
 
 #include <algorithm>
+#include <cmath>
 #include <functional>
 #include <map>
 #include <queue>
@@ -98,6 +100,8 @@ class RooBatchComputeClass : public RooBatchComputeInterface {
    };
 
    void compute(Config const &, Computer computer, std::span<double> output, VarSpan vars, ArgSpan extraArgs) override;
+   void computeExprProgram(Config const &, std::span<const ExprInstr> code, unsigned int stackDepth,
+                           std::span<double> output, VarSpan vars) override;
    double reduceSum(Config const &, InputArr input, size_t n) override;
    ReduceNLLOutput reduceNLL(Config const &, std::span<const double> probas, std::span<const double> weights,
                              std::span<const double> 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 <class F>
+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 <class F>
+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<const ExprInstr> code, unsigned int stackDepth,
+                                              std::span<double> 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<const double> 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<double, double> getLog(double prob, ReduceNLLOutput &out)
 {
    if (prob <= 0.0) {