From 4ca81d81230dd7f9b6100cec0c8648798cef0300 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 18 Oct 2026 09:39:59 +0000 Subject: [PATCH 1/2] [RF] Remove the legacy evaluation backend from the benchmarks The legacy RooFit evaluation backend is deprecated and about to be removed, so the benchmarks stop running and mentioning it: * benchRooFitBackends: drop the *_FitLegacy benchmarks and the unused run configurations. The plot script normalized all timings to the legacy backend; it now shows the speedup of the CUDA backend relative to the CPU backend. * benchRooFitUnbinned: drop the *_FitLegacy benchmarks, including the NumCPU variants. * benchHistFactory: drop the BatchMode_OFF benchmarks. * benchCodeSquashAD: drop the reference minimization with the legacy NLL. The plot script parses the log by fixed line offsets, which shift by one without the legacy NLL creation line. While verifying the new offsets against a log, it turned out that the NLL creation lines of the two codegen variants were assigned the wrong way round; this is fixed as well. * roofitAtlasHiggsBenchmark: drop the legacy NLL slot and the commented-out benchmarks using it. * run_benchmarks.sh, compare_benchmarks.py and the ATLAS benchmark scripts and README no longer run, read or describe a legacy backend. The benchmark names that still contain "BatchMode" are kept unchanged, so that the time series of their results stay comparable. --- root/roofit/atlas-benchmarks/README.md | 6 +-- .../atlasHiggsBackendComparison.C | 2 +- .../atlasHiggsBackendComparison_make_plot.py | 14 +++---- .../roofitAtlasHiggsBenchmark.cxx | 27 +++++-------- .../run_atlasHiggsBackendComparison.sh | 3 -- root/roofit/histfactory/benchHistFactory.cxx | 3 -- .../roofit/RooFitUnBinnedBenchmarks.cxx | 10 ----- root/roofit/roofit/benchCodeSquashAD.cxx | 10 +---- .../roofit/benchCodeSquashAD_make_plot.py | 28 +++++++------- root/roofit/roofit/benchRooFitBackends.cxx | 11 +----- .../roofit/benchRooFitBackends_make_plot.py | 38 +++---------------- root/roofit/roofit/compare_benchmarks.py | 7 +--- root/roofit/roofit/run_benchmarks.sh | 1 - 13 files changed, 43 insertions(+), 117 deletions(-) diff --git a/root/roofit/atlas-benchmarks/README.md b/root/roofit/atlas-benchmarks/README.md index 2b08cd924..604d3e314 100644 --- a/root/roofit/atlas-benchmarks/README.md +++ b/root/roofit/atlas-benchmarks/README.md @@ -33,10 +33,8 @@ Notes: - The `time` builtin has to stay wrapped around `root` itself: its `user` line is the total bar height, and the measured steps are subtracted from it to get the "Other" slice. -- Out of the box this gives two bars, `cpu` and `codegen`. The `legacy` backend - no longer exists in recent ROOT; on an older ROOT, add it to the loop in - `run_atlasHiggsBackendComparison.sh` to get a third bar. The plot script drops - any backend whose log is absent, so no other change is needed. +- Out of the box this gives two bars, `cpu` and `codegen`. The plot script + drops any backend whose log is absent. - The plot script takes the log directory as an argument, so several runs can be compared: `python3 atlasHiggsBackendComparison_make_plot.py atlas some-other-run` writes one PDF per directory, named after it. diff --git a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C index cc6c9dfef..d84bab350 100644 --- a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C +++ b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C @@ -65,7 +65,7 @@ void atlasHiggsBackendComparison(const char *evalBackend = "cpu", bool doHessian *data, GlobalObservables(*globObs), Offset(true), Optimize(2), EvalBackend(evalBackend))}; // Only the codegen NLL is a RooEvaluatorWrapper, so this cast is not valid - // for the legacy and cpu backends. + // for the cpu backend. if (isCodegen && doHessian) { static_cast(*nll).generateHessian(); } diff --git a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py index b39c8d2d2..7a814fb84 100644 --- a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py +++ b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py @@ -18,7 +18,7 @@ nargs="*", default=["atlas"], metavar="DIR", - help="directory holding {legacy,cpu,codegen}.log (default: atlas)", + help="directory holding {cpu,codegen}.log (default: atlas)", ) parser.add_argument( "--json", @@ -122,13 +122,11 @@ def parse_output(filename): # Bars appear left to right in this order. A backend whose log is absent is -# dropped rather than being an error, so a run without the legacy backend simply -# yields a two-bar plot. -backend_names = {"RooFit legacy": "legacy", "RooFit": "cpu", "RooFit AD": "codegen"} - -# The legacy backend no longer exists in recent ROOT, so out of the box this -# produces a two-bar plot. To compare ROOT versions instead of backends, relabel: -# backend_names = {"ROOT 6.32": "legacy", "ROOT 6.40": "cpu", "ROOT 6.40 AD": "codegen"} +# dropped rather than being an error, so a run without one of the backends +# simply yields fewer bars. To compare ROOT versions instead of backends, +# relabel, e.g.: +# backend_names = {"ROOT 6.40": "cpu", "ROOT 6.40 AD": "codegen"} +backend_names = {"RooFit": "cpu", "RooFit AD": "codegen"} # One row per (run, backend, step) for --json: a tidy table stays easy to # concatenate across runs and to pivot downstream. diff --git a/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx b/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx index 41f2c31ae..2d5266fbc 100644 --- a/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx +++ b/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx @@ -39,7 +39,6 @@ struct BenchmarkData { RooAbsPdf *pdf = nullptr; std::unique_ptr ownedPdf; RooAbsData *data = nullptr; - std::unique_ptr nll; std::unique_ptr batchedNll; std::unique_ptr codegenNll; }; @@ -52,9 +51,8 @@ BenchmarkData &bmdata() static void benchCreateNLL(benchmark::State &state) { - const std::string evalBackend = state.range(0) == 2 ? "codegen" : (state.range(0) == 1 ? "cpu" : "legacy"); - auto &nllPtr = - state.range(0) == 2 ? bmdata().codegenNll : (state.range(0) == 1 ? bmdata().batchedNll : bmdata().nll); + const std::string evalBackend = state.range(0) == 2 ? "codegen" : "cpu"; + auto &nllPtr = state.range(0) == 2 ? bmdata().codegenNll : bmdata().batchedNll; auto mc = static_cast(bmdata().ws->obj("ModelConfig")); RooArgSet const *globObs = mc->GetGlobalObservables(); @@ -74,8 +72,7 @@ static void benchCreateNLL(benchmark::State &state) static void benchEvaluateNLL(benchmark::State &state) { bool kickAllParameters = state.range(0); - auto &nllPtr = - state.range(1) == 2 ? bmdata().codegenNll : (state.range(1) == 1 ? bmdata().batchedNll : bmdata().nll); + auto &nllPtr = state.range(1) == 2 ? bmdata().codegenNll : bmdata().batchedNll; RooArgSet parameters; nllPtr->getParameters(nullptr, parameters); @@ -96,8 +93,7 @@ static void benchEvaluateNLL(benchmark::State &state) static void benchMinimizeNLL(benchmark::State &state) { - auto &nllPtr = - state.range(0) == 2 ? bmdata().codegenNll : (state.range(0) == 1 ? bmdata().batchedNll : bmdata().nll); + auto &nllPtr = state.range(0) == 2 ? bmdata().codegenNll : bmdata().batchedNll; RooArgSet parameters; nllPtr->getParameters(nullptr, parameters); @@ -115,28 +111,23 @@ static void benchMinimizeNLL(benchmark::State &state) parameters.assign(initialParams); } -// BENCHMARK(benchCreateNLL)->Name("createNLL")->Args({0})->Unit(kSecond)->Iterations(1); // BENCHMARK(benchCreateNLL)->Name("createNLL_CPU")->Args({1})->Unit(kSecond)->Iterations(1); BENCHMARK(benchCreateNLL)->Name("createNLL_CodeGenAD")->Args({2})->Unit(kSecond)->Iterations(1); -// BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL")->Args({1, 0})->Unit(kMillisecond); // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CPU")->Args({1, 1})->Unit(kMillisecond); // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CodeGenAD")->Args({1, 2})->Unit(kMillisecond); -// BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_SingleKick")->Args({0, 0})->Unit(kMillisecond); // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CPU_SingleKick")->Args({0, 1})->Unit(kMillisecond); // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CodeGenAD_SingleKick")->Args({0, 2})->Unit(kMillisecond); -// BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL")->Args({0})->Unit(kSecond)->Iterations(1); // BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL_CPU")->Args({1})->Unit(kSecond)->Iterations(1); BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL_CodeGenAD")->Args({2})->Unit(kSecond)->Iterations(1); // The channels 221 to 231 inclusive of the full combination workspace are // unfortunately corrupt. They contain RooAddPdfs that are affected by the // notorious server-proxy-desyncing that can happen if the RooFit frameworks do -// the wrong thing. Th new the EvalBackend uses the client-server links to build -// the computation graph for evaluation, and the old RooFit uses the proxies -// (the client-server links are only used to the dirty flag propagation in the -// old RooFit). As the servers and proxies are out of sync in some channels, we -// have no idea what is actually correct and hence these channels are -// unsuitable for validation and benchmarking. They are replaced here by +// the wrong thing. The evaluation backends use the client-server links to +// build the computation graph for evaluation, while the proxies are what the +// pdfs were built with. As the servers and proxies are out of sync in some +// channels, we have no idea what is actually correct and hence these channels +// are unsuitable for validation and benchmarking. They are replaced here by // uniform pdfs. void maskBrokenChannels() { diff --git a/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh b/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh index c7727b33b..3d7add3a6 100755 --- a/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh +++ b/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh @@ -8,9 +8,6 @@ # # The `time` builtin has to wrap root itself, because the plot script reads its # "user" line as the total bar height. -# -# The legacy backend is not run here, as it no longer exists in recent ROOT. -# Add it to the loop below on an older ROOT to get a third bar. mkdir -p atlas diff --git a/root/roofit/histfactory/benchHistFactory.cxx b/root/roofit/histfactory/benchHistFactory.cxx index a9eae01dd..1831d38d5 100644 --- a/root/roofit/histfactory/benchHistFactory.cxx +++ b/root/roofit/histfactory/benchHistFactory.cxx @@ -81,11 +81,8 @@ static void benchHistFactory001(benchmark::State &state) auto const unit = benchmark::kMillisecond; -auto Legacy = static_cast(RooFit::EvalBackend::Value::Legacy); auto Cpu = static_cast(RooFit::EvalBackend::Value::Cpu); -BENCHMARK(benchHistFactory001)->Unit(unit)->Args({1, Legacy})->Name("hf001__BinnedFitOptimization_ON___BatchMode_OFF"); -BENCHMARK(benchHistFactory001)->Unit(unit)->Args({0, Legacy})->Name("hf001__BinnedFitOptimization_OFF__BatchMode_OFF"); BENCHMARK(benchHistFactory001)->Unit(unit)->Args({1, Cpu})->Name("hf001__BinnedFitOptimization_ON___BatchMode_ON_"); BENCHMARK(benchHistFactory001)->Unit(unit)->Args({0, Cpu})->Name("hf001__BinnedFitOptimization_OFF__BatchMode_ON_"); diff --git a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx b/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx index 58f8ecf57..38ee3cb72 100644 --- a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx +++ b/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx @@ -236,31 +236,21 @@ static void BDecayDoubleGauss(benchmark::State &state) int nEvents = 100000; const auto unit = benchmark::kMillisecond; -auto Legacy = static_cast(RooFit::EvalBackend::Value::Legacy); auto Cpu = static_cast(RooFit::EvalBackend::Value::Cpu); auto Cuda = static_cast(RooFit::EvalBackend::Value::Cuda); #define ARGS UseRealTime()->Unit(unit) -BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS; -BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS; -BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS; BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; #ifdef DO_BENCH_ROOFIT_CUDA BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; #endif -BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS; -BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS; -BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS; BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; #ifdef DO_BENCH_ROOFIT_CUDA BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; #endif -BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS; -BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS; -BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS; BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; #ifdef DO_BENCH_ROOFIT_CUDA BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; diff --git a/root/roofit/roofit/benchCodeSquashAD.cxx b/root/roofit/roofit/benchCodeSquashAD.cxx index ebc9648c1..17152fdb7 100644 --- a/root/roofit/roofit/benchCodeSquashAD.cxx +++ b/root/roofit/roofit/benchCodeSquashAD.cxx @@ -63,7 +63,6 @@ std::unique_ptr createModel(RooRealVar &x, std::string const &channel } enum backend { - Reference, BatchMode, CodeSquashNumDiff, CodeSquashAD @@ -126,7 +125,6 @@ static void BM_RooFuncWrapper_ManyParams_Minimization(benchmark::State &state) std::cout << "nparams: " << params.size() << std::endl; - std::unique_ptr nllRef{model.createNLL(data, EvalBackend::Legacy(), Offset("initial"))}; std::unique_ptr nllRefBatch{model.createNLL(data, EvalBackend::Cpu(), Offset("initial"))}; std::unique_ptr nllFunc{model.createNLL(data, EvalBackend::Codegen(), Offset("initial"))}; std::unique_ptr nllFuncNoGrad{model.createNLL(data, EvalBackend::CodegenNoGrad(), Offset("initial"))}; @@ -134,9 +132,7 @@ static void BM_RooFuncWrapper_ManyParams_Minimization(benchmark::State &state) std::unique_ptr m = nullptr; int code = state.range(0); - if (code == backend::Reference) { - m = std::make_unique(*nllRef); - } else if (code == backend::CodeSquashNumDiff) { + if (code == backend::CodeSquashNumDiff) { m = std::make_unique(*nllFuncNoGrad); } else if (code == backend::BatchMode) { m = std::make_unique(*nllRefBatch); @@ -172,10 +168,6 @@ int main(int argc, char **argv) range.insert(range.end(), rangeHigh.begin(), rangeHigh.end()); #endif - // Run the minimization with the reference NLL - RooFitADBenchmarksUtils::doBenchmarks("NllReferenceMinimization", Reference, - BM_RooFuncWrapper_ManyParams_Minimization, range, 1); - // Run the minimization with the reference NLL (BatchMode) RooFitADBenchmarksUtils::doBenchmarks("NllBatchModeMinimization", BatchMode, BM_RooFuncWrapper_ManyParams_Minimization, range, 1); diff --git a/root/roofit/roofit/benchCodeSquashAD_make_plot.py b/root/roofit/roofit/benchCodeSquashAD_make_plot.py index 7ed3f33d8..4979c6221 100644 --- a/root/roofit/roofit/benchCodeSquashAD_make_plot.py +++ b/root/roofit/roofit/benchCodeSquashAD_make_plot.py @@ -56,20 +56,22 @@ def is_relevant(line): lines = list(filter(is_relevant, lines)) -names = ["NllReferenceMinimization", "NllBatchModeMinimization", "NllCodeSquash_NumDiff", "NllCodeSquash_AD"] +names = ["NllBatchModeMinimization", "NllCodeSquash_NumDiff", "NllCodeSquash_AD"] backends = { - "NllReferenceMinimization": "legacy", "NllBatchModeMinimization": "RooFit", "NllCodeSquash_NumDiff": "Hardcoded", "NllCodeSquash_AD": "RooFit AD", } +# Index of the "Creation of NLL object took" line of each backend in the block +# of relevant lines per benchmark (see analyze_block). The NLL objects are +# created in the order cpu, codegen (followed by its three jitting lines), +# codegen_no_grad. nll_line = { - "NllReferenceMinimization": 1, - "NllBatchModeMinimization": 2, - "NllCodeSquash_NumDiff": 3, - "NllCodeSquash_AD": 7, + "NllBatchModeMinimization": 1, + "NllCodeSquash_AD": 2, + "NllCodeSquash_NumDiff": 6, } @@ -77,7 +79,7 @@ def analyze_block(lines): # print(lines) data = dict() - name = lines[12].split("/")[0] + name = lines[11].split("/")[0] data["backend"] = backends[name] @@ -86,17 +88,17 @@ def line_to_seconds(l): return time_string_to_seconds(" ".join(l.split(" ")[-2:])) data["create_nll"] = line_to_seconds(lines[nll_line[name]]) - data["seeding"] = line_to_seconds(lines[9]) - data["migrad"] = line_to_seconds(lines[10]) + data["seeding"] = line_to_seconds(lines[8]) + data["migrad"] = line_to_seconds(lines[9]) data["nparams"] = int(lines[0].split(" ")[-1]) - 1 - data["fval"] = float(lines[11].split(" ")[-1]) + data["fval"] = float(lines[10].split(" ")[-1]) - data["jitting"] = sum(line_to_seconds(lines[i]) for i in [4, 5, 6]) + data["jitting"] = sum(line_to_seconds(lines[i]) for i in [3, 4, 5]) return data -m_block = 13 +m_block = 12 datas = [] @@ -107,7 +109,7 @@ def line_to_seconds(l): df = pd.DataFrame(datas) -df = df.query("nparams > 10 and nparams < 2500 and backend != 'legacy'") +df = df.query("nparams > 10 and nparams < 2500") dfs = dict() diff --git a/root/roofit/roofit/benchRooFitBackends.cxx b/root/roofit/roofit/benchRooFitBackends.cxx index 887d788a7..91836c316 100644 --- a/root/roofit/roofit/benchRooFitBackends.cxx +++ b/root/roofit/roofit/benchRooFitBackends.cxx @@ -25,7 +25,7 @@ int printLevel = 0; size_t nEvents = 10000; const auto minimizerName = "Minuit2"; -enum RunConfig_t { runScalar, runCpu, fitScalar, fitCpu, fitCuda, fitCodegen, fitCodegenNoGrad }; +enum RunConfig_t { fitCpu, fitCuda, fitCodegen, fitCodegenNoGrad }; void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data) { @@ -38,9 +38,7 @@ void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data) using namespace RooFit; for (auto _ : state) { - if (runConfig == fitScalar) { - pdf.fitTo(data, EvalBackend::Legacy(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1)); - } else if (runConfig == fitCpu) { + if (runConfig == fitCpu) { pdf.fitTo(data, EvalBackend::Cpu(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1)); } else if (runConfig == fitCuda) { pdf.fitTo(data, EvalBackend::Cuda(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1)); @@ -326,7 +324,6 @@ auto const unit = benchmark::kMillisecond; #define CUDA_ONLY(code) #endif -BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitLegacy")->Args({fitScalar}); BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCPU")->Args({fitCpu}); BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegen")->Args({fitCodegen}); BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegenNoGrad")->Args({fitCodegenNoGrad}); @@ -335,17 +332,14 @@ CUDA_ONLY(BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCUDA")->Args({fitC // No codegen variants for this model: normalizing the Gaussian over both x // and sigma requires a multi-dimensional numeric integral, which the codegen // backend doesn't support. -BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitLegacy")->Args({fitScalar}); BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCPU")->Args({fitCpu}); CUDA_ONLY(BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCUDA")->Args({fitCuda})); -BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitLegacy")->Args({fitScalar}); BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCPU")->Args({fitCpu}); BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegen")->Args({fitCodegen}); BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegenNoGrad")->Args({fitCodegenNoGrad}); CUDA_ONLY(BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCUDA")->Args({fitCuda})); -BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitLegacy")->Unit(benchmark::kMillisecond)->Args({fitScalar}); BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu}); BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen}); BENCHMARK(benchProdPdf) @@ -358,7 +352,6 @@ CUDA_ONLY(BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCUDA")->Unit(ben // Watch out with the result from these benchmarks: if there are evaluation // errors during the minimization, the time differences will mostly come from // how the backends are handling those. This might or might not be wanted. -BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitLegacy")->Unit(benchmark::kMillisecond)->Args({fitScalar}); BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu}); BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen}); BENCHMARK(benchModel) diff --git a/root/roofit/roofit/benchRooFitBackends_make_plot.py b/root/roofit/roofit/benchRooFitBackends_make_plot.py index 0480b37d8..65ab1f2ed 100644 --- a/root/roofit/roofit/benchRooFitBackends_make_plot.py +++ b/root/roofit/roofit/benchRooFitBackends_make_plot.py @@ -34,8 +34,8 @@ # "FitModel", # ] -backends = ["CPU", "CUDA"] -reference_backend = "Legacy" +backends = ["CUDA"] +reference_backend = "CPU" timings = {} @@ -87,7 +87,7 @@ ROOT.gStyle.SetHistMinimumZero() h1b = ROOT.TH1F( - "h1b", "RooFit: speedup in benchmark fits with BatchMode() relative to old RooFit (1 million events)", nx, 0, nx + "h1b", "RooFit: speedup in benchmark fits with the CUDA backend relative to the CPU backend (1 million events)", nx, 0, nx ) h1b.SetFillColor(9) h1b.SetBarWidth(0.8 / nb) @@ -102,41 +102,13 @@ # h1b.GetXaxis().SetLabelOffset(0.1) for i in range(nx): - h1b.SetBinContent(i + 1, timings["CPU"][i]) + h1b.SetBinContent(i + 1, timings["CUDA"][i]) h1b.GetXaxis().SetBinLabel(i + 1, labels[i]) h1b.Draw("b") -h2b = ROOT.TH1F("h2b", "h2b", nx, 0, nx) - -h2b.SetFillColor(46) -h2b.SetBarWidth(0.8 / nb) -h2b.SetBarOffset(0.1 + 2 * 0.8 / nb) -h2b.SetStats(0) -h2b.SetMinimum(-1) -h2b.SetMaximum(y_max) -for i in range(nx): - h2b.SetBinContent(i + 1, 0.0) -h2b.Draw("b same") - -h3b = ROOT.TH1F("h3b", "h3b", nx, 0, nx) - -h3b.SetFillColor(8) -h3b.SetBarWidth(0.8 / nb) -h3b.SetBarOffset(0.1 + 1 * 0.8 / nb) -h3b.SetStats(0) -h3b.SetMinimum(-1) -h3b.SetMaximum(y_max) -for i in range(nx): - h3b.SetBinContent(i + 1, timings["CUDA"][i]) - -h3b.Draw("b same") - legend = ROOT.TLegend(0.1, 0.7, 0.46, 0.9) -# legend->SetHeader("The Legend Title","C"); // option "C" allows to center the heade -legend.AddEntry(h1b, "CPU (AMD Ryzen 9 3900)", "f") -# legend.AddEntry(h2b, "Vectorized (24 threads)", "f") -legend.AddEntry(h3b, "CUDA (RTX 3070 gaming GPU)", "f") +legend.AddEntry(h1b, "CUDA relative to CPU (single thread)", "f") legend.Draw() cb.SaveAs(base_name + ".png") diff --git a/root/roofit/roofit/compare_benchmarks.py b/root/roofit/roofit/compare_benchmarks.py index d92734201..5577880d7 100644 --- a/root/roofit/roofit/compare_benchmarks.py +++ b/root/roofit/roofit/compare_benchmarks.py @@ -14,7 +14,6 @@ def read_csv(csv_file_path): codegen_df = read_csv("out_codegen.csv") codegen_nograd_df = read_csv("out_codegen_ngrad.csv") -legacy_df = read_csv("out_legacy.csv") cpu_df = read_csv("out_cpu.csv") # Plotting @@ -28,17 +27,15 @@ def read_csv(csv_file_path): codegen_time = codegen_df.loc[codegen_df["name"] == benchmark, "real_time"] codegen_nograd_time = codegen_nograd_df.loc[codegen_nograd_df["name"] == benchmark, "real_time"] cpu_time = cpu_df.loc[cpu_df["name"] == benchmark, "real_time"] - legacy_time = legacy_df.loc[legacy_df["name"] == benchmark, "real_time"] plt.bar(x[i] - 0.10, codegen_time, width=0.15, align="center", label="codegen", color="lightblue") plt.bar(x[i], codegen_nograd_time, width=0.15, align="edge", label="codegen_nograd", color="navy") plt.bar(x[i] + 0.15, cpu_time, width=0.15, align="edge", label="cpu", color="cyan") - plt.bar(x[i] + 0.30, legacy_time, width=0.15, align="edge", label="legacy", color="gray") # Customize legend -legend_labels = ["codegen", "codegen_nograd", "cpu", "legacy"] -legend_colors = ["lightblue", "navy", "cyan", "gray"] +legend_labels = ["codegen", "codegen_nograd", "cpu"] +legend_colors = ["lightblue", "navy", "cyan"] legend_handles = [plt.Rectangle((0, 0), 1, 1, color=color) for color in legend_colors] plt.legend(legend_handles, legend_labels) diff --git a/root/roofit/roofit/run_benchmarks.sh b/root/roofit/roofit/run_benchmarks.sh index 308388f86..b5b9da1ff 100755 --- a/root/roofit/roofit/run_benchmarks.sh +++ b/root/roofit/roofit/run_benchmarks.sh @@ -20,7 +20,6 @@ run_benchmark() { # Run benchmarks run_benchmark "./benchRooFitBinned -b codegen --benchmark_out_format=csv --benchmark_out=out_codegen.csv" "out_codegen.csv" run_benchmark "./benchRooFitBinned -b codegen_no_grad --benchmark_out_format=csv --benchmark_out=out_codegen_ngrad.csv" "out_codegen_ngrad.csv" -run_benchmark "./benchRooFitBinned -b legacy --benchmark_out_format=csv --benchmark_out=out_legacy.csv" "out_legacy.csv" run_benchmark "./benchRooFitBinned -b cpu --benchmark_out_format=csv --benchmark_out=out_cpu.csv" "out_cpu.csv" # Run Python script From 0f8b53420aeb765e9be006790825e6bd38f90204 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 18 Oct 2026 13:16:50 +0000 Subject: [PATCH 2/2] [RF] Consolidate the unbinned fit backend benchmarks Merge benchRooFitUnbinned into benchRooFitBackends, which now covers all unbinned models with one protocol: the likelihood is created outside of the timed loop, and only Migrad with strategy 0 is timed, after one untimed warm-up fit and with a fresh minimizer for every fit. * Models from both executables in one file, 100k events by default (--events=N option), fixed random seed for the datasets, so that results are reproducible. * Benchmarks for the single-threaded CPU backend, the multi-threaded CPU backend with all threads of the implicit MT pool, CUDA, and codegen where it supports the model and ROOT has clad. * Wall-clock time instead of CPU time, since the multi-threaded entry would otherwise be reported as the sum over all threads. * Counters for the number of likelihood evaluations, the Migrad status and the minimum NLL, to check that the backends do the same work. * The run information (CPU, GPU, ROOT version, RooBatchCompute library, VDT, events, threads) goes into the benchmark context of the JSON output. * The plot script shows the speedup of the multi-threaded CPU and the CUDA backend relative to the single-threaded CPU backend for all models, as horizontal bars with a fixed axis range, together with the run information. It uses the median when the benchmarks were run with --benchmark_repetitions. --- root/roofit/roofit/CMakeLists.txt | 6 - .../roofit/RooFitUnBinnedBenchmarks.cxx | 261 ------- root/roofit/roofit/benchRooFitBackends.cxx | 713 +++++++++++------- .../roofit/benchRooFitBackends_make_plot.py | 293 ++++--- 4 files changed, 626 insertions(+), 647 deletions(-) delete mode 100644 root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx diff --git a/root/roofit/roofit/CMakeLists.txt b/root/roofit/roofit/CMakeLists.txt index 3341bd26e..bc1702ba9 100644 --- a/root/roofit/roofit/CMakeLists.txt +++ b/root/roofit/roofit/CMakeLists.txt @@ -3,18 +3,12 @@ RB_ADD_GBENCHMARK(benchRooFitBinned LABEL long LIBRARIES Core Hist MathCore RIO RooFit RooStats RooFitCore HistFactory) -RB_ADD_GBENCHMARK(benchRooFitUnbinned - RooFitUnBinnedBenchmarks.cxx - LABEL long - LIBRARIES Core Hist MathCore RIO RooFit RooStats RooFitCore HistFactory) - RB_ADD_GBENCHMARK(benchRooFitBackends benchRooFitBackends.cxx LABEL long LIBRARIES Core MathCore RooFitCore RooFit) if(cuda) - target_compile_definitions(benchRooFitUnbinned PRIVATE DO_BENCH_ROOFIT_CUDA) target_compile_definitions(benchRooFitBackends PRIVATE DO_BENCH_ROOFIT_CUDA) endif() diff --git a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx b/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx deleted file mode 100644 index 38ee3cb72..000000000 --- a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx +++ /dev/null @@ -1,261 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -namespace { - -void shiftParameters(const RooArgSet ¶meters) -{ - // We need to change the parameters so we don't start already at the minimum - // (the parameter values for which the dataset was created). This time, we - // don't randomize them because the fit has many parameter and is not stable - // from many initial points. That's why they are only shifted away a little - // bit from their true value. - for (auto *param : static_range_cast(parameters)) { - double val = param->getVal(); - val = val + 0.05 * (param->getMax() - val); - param->setVal(val); - } -} - -std::unique_ptr makeChangeMsgLvl() -{ - auto &msg = RooMsgService::instance(); - msg.reset(); - auto oldKillBelow = msg.globalKillBelow(); - auto out = std::make_unique(RooFit::WARNING); - msg.setGlobalKillBelow(oldKillBelow); - // Keep the FastEvaluations info so we see which classes still need to be supported - msg.addStream(RooFit::INFO, Topic(RooFit::FastEvaluations)); - return out; -} - -} // namespace - -static void BDecayWithMixing(benchmark::State &state) -{ - using namespace RooFit; - - auto changeMsgLevel = makeChangeMsgLvl(); - - gErrorIgnoreLevel = kInfo; - int events = state.range(0); - auto evalBackend = static_cast(state.range(1)); - int cpu = state.range(2); - - RooRealVar dt("dt", "dt", -10, 10); - dt.setBins(40); - // Parameters - RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9); - RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9); - RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0); - RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0); - - RooCategory mixState("mixState", "B0/B0bar mixing state"); - mixState.defineType("mixed", -1); - mixState.defineType("unmixed", 1); - - RooCategory tagFlav("tagFlav", "Flavour of the tagged B0"); - tagFlav.defineType("B0", 1); - tagFlav.defineType("B0bar", -1); - - // Use delta function resolution model - RooTruthModel tm("tm", "truth model", dt); - // Construct Bdecay with mixing - RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided); - // Generate Some Data - std::unique_ptr data{bmix.generate({dt, mixState, tagFlav}, events)}; - - // Save initial parameter snapshot for resetting - RooArgSet params; - RooArgSet paramsInitial; - bmix.getParameters(data->get(), params); - shiftParameters(params); - params.snapshot(paramsInitial); - - // Create NLL - std::unique_ptr nll{bmix.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))}; - RooMinimizer m(*nll); - m.setPrintLevel(-1); - m.setStrategy(0); - m.setLogFile("benchmigradnchanellog"); - - while (state.KeepRunning()) { - m.migrad(); - params.assign(paramsInitial); - } -} - -static void BDecayGaussResolution(benchmark::State &state) -{ - using namespace RooFit; - - auto changeMsgLevel = makeChangeMsgLvl(); - - gErrorIgnoreLevel = kInfo; - int events = state.range(0); - auto evalBackend = static_cast(state.range(1)); - int cpu = state.range(2); - - RooRealVar dt("dt", "dt", -10, 10); - dt.setBins(40); - // Parameters - RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9); - RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9); - RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0); - RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0); - - RooCategory mixState("mixState", "B0/B0bar mixing state"); - mixState.defineType("mixed", -1); - mixState.defineType("unmixed", 1); - - RooCategory tagFlav("tagFlav", "Flavour of the tagged B0"); - tagFlav.defineType("B0", 1); - tagFlav.defineType("B0bar", -1); - - // Use delta function resolution model - RooTruthModel tm("tm", "truth model", dt); - // Construct Bdecay with mixing - RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided); - // Build a gaussian resolution model - RooRealVar bias1("bias1", "bias1", 0, -1., 1.); - RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9); - RooGaussModel gm1("gm1", "gauss model 1", dt, bias1, sigma1); - - // Construct decay(t) (x) gauss1(t) - RooDecay decay_gm1("decay_gm1", "decay", dt, tau, gm1, RooDecay::DoubleSided); - - // Generate Some Data - std::unique_ptr data{decay_gm1.generate({dt, mixState, tagFlav}, events)}; - - // Save initial parameter snapshot for resetting - RooArgSet params; - RooArgSet paramsInitial; - decay_gm1.getParameters(data->get(), params); - shiftParameters(params); - params.snapshot(paramsInitial); - - // Create NLL - std::unique_ptr nll{decay_gm1.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))}; - RooMinimizer m(*nll); - m.setPrintLevel(-1); - m.setStrategy(0); - m.setLogFile("benchmigradnchanellog"); - while (state.KeepRunning()) { - m.migrad(); - params.assign(paramsInitial); - } -} - -static void BDecayDoubleGauss(benchmark::State &state) -{ - using namespace RooFit; - - RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); - - gErrorIgnoreLevel = kInfo; - int events = state.range(0); - auto evalBackend = static_cast(state.range(1)); - int cpu = state.range(2); - - RooRealVar dt("dt", "dt", -10, 10); - dt.setBins(40); - // Parameters - RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9); - RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9); - RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0); - RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0); - - RooCategory mixState("mixState", "B0/B0bar mixing state"); - mixState.defineType("mixed", -1); - mixState.defineType("unmixed", 1); - - RooCategory tagFlav("tagFlav", "Flavour of the tagged B0"); - tagFlav.defineType("B0", 1); - tagFlav.defineType("B0bar", -1); - - // Use delta function resolution model - RooTruthModel tm("tm", "truth model", dt); - // Construct Bdecay with mixing - RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided); - // Build a gaussian resolution model - RooRealVar bias1("bias1", "bias1", 0, -1., 1.); - RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9); - RooGaussModel gm1("gm1", "gauss model 1", dt, bias1, sigma1); - - // Construct decay(t) (x) gauss1(t) - RooDecay decay_gm1("decay_gm1", "decay", dt, tau, gm1, RooDecay::DoubleSided); - - // Build another gaussian resolution model - RooRealVar bias2("bias2", "bias2", 0, -1., 1.); - RooRealVar sigma2("sigma2", "sigma2", 5, 0.5, 10); - RooGaussModel gm2("gm2", "gauss model 2", dt, bias2, sigma2); - // Build a composite resolution model f*gm1+(1-f)*gm2 - RooRealVar gm1frac("gm1frac", "fraction of gm1", 0.5, 0.0, 1.0); - RooAddModel gmsum("gmsum", "sum of gm1 and gm2", {gm1, gm2}, gm1frac); - // Construct decay(t) (x) (f*gm1 + (1-f)*gm2) - RooDecay decay_gmsum("decay_gmsum", "decay", dt, tau, gmsum, RooDecay::DoubleSided); - - // Generate Some Data - std::unique_ptr data{decay_gmsum.generate({dt, mixState, tagFlav}, events)}; - - // Save initial parameter snapshot for resetting - RooArgSet params; - RooArgSet paramsInitial; - decay_gmsum.getParameters(data->get(), params); - shiftParameters(params); - params.snapshot(paramsInitial); - - // Create NLL - std::unique_ptr nll{decay_gmsum.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))}; - RooMinimizer m(*nll); - m.setPrintLevel(-1); - m.setStrategy(0); - m.setLogFile("benchmigradnchanellog"); - while (state.KeepRunning()) { - m.migrad(); - params.assign(paramsInitial); - } -} - -int nEvents = 100000; -const auto unit = benchmark::kMillisecond; - -auto Cpu = static_cast(RooFit::EvalBackend::Value::Cpu); -auto Cuda = static_cast(RooFit::EvalBackend::Value::Cuda); - -#define ARGS UseRealTime()->Unit(unit) - -BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; -#ifdef DO_BENCH_ROOFIT_CUDA -BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; -#endif - -BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; -#ifdef DO_BENCH_ROOFIT_CUDA -BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; -#endif - -BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCPU")->Args({nEvents, Cpu, 1})->ARGS; -#ifdef DO_BENCH_ROOFIT_CUDA -BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS; -#endif - -#undef ARGS - -BENCHMARK_MAIN(); diff --git a/root/roofit/roofit/benchRooFitBackends.cxx b/root/roofit/roofit/benchRooFitBackends.cxx index 91836c316..eaa90481f 100644 --- a/root/roofit/roofit/benchRooFitBackends.cxx +++ b/root/roofit/roofit/benchRooFitBackends.cxx @@ -10,376 +10,513 @@ * listed in LICENSE (http://roofit.sourceforge.net/license.txt) */ +// Benchmarks of unbinned maximum likelihood fits with the different RooFit +// evaluation backends. Each benchmark builds a model, generates a dataset of +// nEvents events from it, moves the parameters away from the generation +// values, and then times Migrad with strategy 0 on a likelihood that was +// created outside of the timed loop. The benchmarks are named +// "_Fit", which the script benchRooFitBackends_make_plot.py +// relies on. The CPU backend is benchmarked single-threaded (CPU) and with +// the multi-threaded evaluation on all threads of ROOT's implicit +// multi-threading pool (CPUMT). The benchmarks report wall-clock time, which +// is what matters when comparing thread counts. +// +// The information that is needed to label the results (CPU and GPU, ROOT +// version, the RooBatchCompute library that was loaded, VDT, number of +// events and threads) is written to the context of the benchmark output, +// which the JSON output format includes. +// +// The number of events per fit is 100k by default and can be changed with the +// "--events=N" option. Each benchmark does one untimed warm-up fit before the +// timed iterations. The slower fits take only one iteration, so run with +// "--benchmark_repetitions=3" to get a spread estimate; the plot script then +// uses the median. The number of likelihood evaluations, the Migrad status +// and the minimum NLL are reported as counters to check that the backends do +// the same work. + +#include #include +#include +#include #include +#include #include -#include +#include #include +#include +#include +#include +#include #include #include +#include #include +#include +#include +#include + #include -int printLevel = 0; -size_t nEvents = 10000; -const auto minimizerName = "Minuit2"; +#include +#include +#include +#include +#include +#include -enum RunConfig_t { fitCpu, fitCuda, fitCodegen, fitCodegenNoGrad }; +namespace { -void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data) -{ - RunConfig_t runConfig = static_cast(state.range(0)); - - RooArgSet params; - pdf.getParameters(data.get(), params); - RooArgSet paramsInitial; - params.snapshot(paramsInitial); +std::size_t nEvents = 100000; // can be set with the --events option - using namespace RooFit; - for (auto _ : state) { - if (runConfig == fitCpu) { - pdf.fitTo(data, EvalBackend::Cpu(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1)); - } else if (runConfig == fitCuda) { - pdf.fitTo(data, EvalBackend::Cuda(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1)); - } else if (runConfig == fitCodegen) { - pdf.fitTo(data, EvalBackend::Codegen(), Minimizer(minimizerName), PrintLevel(printLevel - 1), - PrintEvalErrors(-1)); - } else if (runConfig == fitCodegenNoGrad) { - pdf.fitTo(data, EvalBackend::CodegenNoGrad(), Minimizer(minimizerName), PrintLevel(printLevel - 1), - PrintEvalErrors(-1)); - } - state.PauseTiming(); - params.assign(paramsInitial); - state.ResumeTiming(); - } +/// Silence RooFit, but keep the FastEvaluations info so we see which classes +/// still need to be supported by the evaluation backends. +std::unique_ptr makeChangeMsgLvl() +{ + auto &msg = RooMsgService::instance(); + msg.reset(); + auto oldKillBelow = msg.globalKillBelow(); + auto out = std::make_unique(RooFit::WARNING); + msg.setGlobalKillBelow(oldKillBelow); + msg.addStream(RooFit::INFO, RooFit::Topic(RooFit::FastEvaluations)); + return out; } -static void benchProdPdf(benchmark::State &state) +/// Move the parameters away from the values that the dataset was generated +/// with, so that the fit doesn't start at the minimum. The parameters are not +/// randomized, because fits with many parameters are not stable from arbitrary +/// starting points: they are only shifted towards their upper limits by the +/// given fraction of the distance. +void shiftParameters(RooArgSet const ¶meters, double fraction) { - RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING); - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range - - RooWorkspace w; - w.factory("Gaussian::g1(x[0,20],mx1[10,0,20],sx1[3,0.1,10])"); - w.factory("Gaussian::g2(x,mx2[10,0,20],sx2[1,0.1,10])"); - w.factory("Gaussian::g3(y[0,20],my[10,0,20],sy[2,0.1,10])"); - w.factory("Polynomial::pol(y,a[-0.01,-0.05,0.1])"); - w.factory("Gamma::gpdf(z[0,10],g[2,0.1,10],b[2,0.1,10],m0[-1,-20,0])"); - - w.factory("SUM::pdfx(fx[0.2,0,1]*g1,g2)"); - w.factory("SUM::pdfy(fy[0.4,0.,1.]*g3,pol)"); - w.factory("PROD::model(pdfx,pdfy,gpdf)"); - - // set coeffients - // w.var("ns1")->setVal(0.2*nEvents); - // w.var("ns2")->setVal(0.3*nEvents); - // w.var("ng2")->setVal(0.1*nEvents); - // w.var("ng3")->setVal(0.1*nEvents); - // w.var("npol")->setVal(0.3*nEvents); - - auto &pdf = *(w.pdf("model")); - auto &x = *(w.var("x")); - auto &y = *(w.var("y")); - auto &z = *(w.var("z")); - - std::unique_ptr data{pdf.generate({x, y, z}, nEvents)}; - - RooArgSet parameters; - pdf.getParameters(data->get(), parameters); - // We need to change the parameters so we don't start already at the minimum - // (the parameter values for which the dataset was created). This time, we - // don't randomize them because the fit has many parameter and is not stable - // from many initial points. That's why they are only shifted away a little - // bit from their true value. for (auto *param : static_range_cast(parameters)) { double val = param->getVal(); - val = val + 0.05 * (param->getMax() - val); - param->setVal(val); + param->setVal(val + fraction * (param->getMax() - val)); } - - runFitBenchmark(state, pdf, *data); } -class GausModel { - -private: - std::unique_ptr _w; - std::unique_ptr _data; - RooAbsPdf *_pdf; - std::string _obsName; - std::vector _results; - -public: - GausModel(std::string const &obsName, size_t nEvts) - { - - RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING); - - _w = std::make_unique("w"); - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range - RooRealVar x("x", "x", -1.5, 40.5); - - RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.); - RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.); +/// Generate the dataset for a benchmark. The random generator is seeded before +/// each generation, such that every benchmark fits the same dataset no matter +/// in which order or how often the benchmarks run. Otherwise, the +/// minimization path and hence the timing would vary from run to run, +/// especially for the models that produce evaluation errors during the fit. +std::unique_ptr generateData(RooAbsPdf &pdf, RooArgSet const &observables) +{ + RooRandom::randomGenerator()->SetSeed(1337); + return std::unique_ptr{pdf.generate(observables, nEvents)}; +} - // Build gaussian p.d.f in terms of x,mean and sigma - RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma); +/// Set the parameters to uniformly distributed random values in their ranges, +/// for models with few parameters where the fit converges from anywhere. +void randomizeParameters(RooArgSet const ¶meters, unsigned long seed) +{ + auto random = RooRandom::randomGenerator(); + random->SetSeed(seed); + for (auto *param : static_range_cast(parameters)) { + param->setVal(param->getMin() + random->Uniform() * (param->getMax() - param->getMin())); + } +} - _w->import(gauss); +/// Create the likelihood with the backend and the number of threads requested +/// by the benchmark arguments and time Migrad with strategy 0, starting from +/// the same parameter values (and errors, i.e. Minuit step sizes) in each +/// iteration. The parameters are expected to be moved away from the minimum +/// by the caller. A thread count of zero means all threads of ROOT's implicit +/// multi-threading pool, which main() enables. +void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data) +{ + using namespace RooFit; - _w->defineSet("obs", obsName.c_str()); + auto evalBackend = static_cast(state.range(0)); + auto nThreads = static_cast(state.range(1)); + if (nThreads == 0) { + nThreads = ROOT::GetThreadPoolSize(); + } + state.counters["num_threads"] = nThreads; - _pdf = _w->pdf("gauss"); + RooArgSet params; + pdf.getParameters(data.get(), params); + RooArgSet paramsInitial; + params.snapshot(paramsInitial); - generateData(nEvts); + std::unique_ptr nll{pdf.createNLL(data, EvalBackend(evalBackend), NumCPU(nThreads))}; + + // Every fit gets a fresh minimizer: the evaluation error recovery state of + // the minimizer function (maximum FCN value, number of invalid + // evaluations, function offset) persists across Migrad calls, so a second + // fit with the same minimizer can take a different path than the first one + // for models with evaluation errors. Creating the minimizer is not timed. + auto makeMinimizer = [&]() { + auto m = std::make_unique(*nll); + m->setPrintLevel(-1); + m->setPrintEvalErrors(-1); + m->setStrategy(0); + return m; + }; + + // One untimed fit first: the first likelihood evaluations pay for one-time + // work like loading libraries or plugins, spinning up the thread pool, or + // loading the CUDA kernels. Google Benchmark reports the first iteration + // as the result if it is slow enough, so without the warm-up that cost + // would end up in the timing of the slow fits only. + makeMinimizer()->migrad(); + params.assign(paramsInitial); + + auto m = makeMinimizer(); + for (auto _ : state) { + const int status = m->migrad(); + state.PauseTiming(); + // The counters show whether the backends do the same work and whether + // the fit succeeded. They are the same in every iteration, so the last + // one is reported. + state.counters["nll_evals"] = m->evalCounter(); + state.counters["migrad_status"] = status; + state.counters["min_nll"] = nll->getVal(); + params.assign(paramsInitial); + m = makeMinimizer(); + state.ResumeTiming(); } - RooAbsPdf &pdf() { return *_pdf; } - RooAbsData &data() { return *_data; } +} - void randomiseParameters(unsigned long seed = 0) +// Common objects of the B decay models. +struct BDecayVars { + RooRealVar dt{"dt", "dt", -10, 10}; + RooRealVar dm{"dm", "delta m(B0)", 0.472, 0.1, 0.9}; + RooRealVar tau{"tau", "tau (B0)", 1.547, 1.3, 1.9}; + RooRealVar w{"w", "flavour mistag rate", 0.1, 0.0, 1.0}; + RooRealVar dw{"dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0}; + RooCategory mixState{"mixState", "B0/B0bar mixing state"}; + RooCategory tagFlav{"tagFlav", "Flavour of the tagged B0"}; + + BDecayVars() { - RooArgSet parameters; - _pdf->getParameters(_data->get(), parameters); - auto random = RooRandom::randomGenerator(); - if (seed != 0) - random->SetSeed(seed); - - for (auto param : parameters) { - auto par = static_cast(param); - const double uni = random->Uniform(); - const double min = par->getMin(); - const double max = par->getMax(); - par->setVal(min + uni * (max - min)); - } + dt.setBins(40); + mixState.defineType("mixed", -1); + mixState.defineType("unmixed", 1); + tagFlav.defineType("B0", 1); + tagFlav.defineType("B0bar", -1); } - void generateData(size_t nEvts) + void generateAndRun(benchmark::State &state, RooAbsPdf &pdf) { - _data = std::unique_ptr(_pdf->generate(*_w->set("obs"), nEvts)); - // allocate here output vector - _results = std::vector(nEvts); + std::unique_ptr data{generateData(pdf, {dt, mixState, tagFlav})}; + RooArgSet params; + pdf.getParameters(data->get(), params); + shiftParameters(params, 0.05); + runFitBenchmark(state, pdf, *data); } }; +} // namespace + +/// Gaussian with one observable. static void benchFitGauss(benchmark::State &state) { - GausModel model("x", nEvents); - // We need to randomize the parameters so we don't start already at the - // minimum (the parameter values for which the dataset was created). - model.randomiseParameters(1337); - runFitBenchmark(state, model.pdf(), model.data()); + auto changeMsgLevel = makeChangeMsgLvl(); + + RooRealVar x("x", "x", -1.5, 40.5); + RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.); + RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.); + RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma); + + std::unique_ptr data{generateData(gauss, {x})}; + randomizeParameters({mean, sigma}, 1337); + runFitBenchmark(state, gauss, *data); } +/// Gaussian with both x and sigma as observables, i.e. a conditional pdf. static void benchFitGaussXSigma(benchmark::State &state) { - GausModel model("x,sigma", nEvents); - // We need to randomize the parameters so we don't start already at the - // minimum (the parameter values for which the dataset was created). - model.randomiseParameters(1337); - runFitBenchmark(state, model.pdf(), model.data()); -} - -class TestModelPdf { - -private: - std::unique_ptr _ws; - std::unique_ptr _data; - RooAbsPdf *_pdf; - std::string _obsName; - std::vector _results; - -public: - TestModelPdf(size_t nEvts) - { - - if (printLevel == 0) - RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING); - - _ws = std::make_unique("w"); - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range - RooRealVar x("x", "x", -1.5, 40.5); - - RooRealVar mean("mean", "mean of gaussian", 20., -10, 30); - RooRealVar sigma("sigma", "width of gaussian", 4., 0.5, 10); - - // Build gaussian p.d.f in terms of x,mean and sigma - RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma); - - RooRealVar c1("c1", "Decay constant", -0.5, -10, -0.001); - RooExponential ex{"Expo", "Exponential PDF", x, c1}; - - RooRealVar fractionGaus("fractionGaus", "Fraction of Gauss component", 0.5, 0., 1.); - RooAddPdf addPdf("SumGausExpo", "Sum of Gaus and Exponential", {gauss, ex}, fractionGaus); - // to avoid a warning when computing the unnormalized RooAddPdf values - addPdf.fixCoefNormalization(x); + auto changeMsgLevel = makeChangeMsgLvl(); - _ws->import(addPdf); + RooRealVar x("x", "x", -1.5, 40.5); + RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.); + RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.); + RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma); - _ws->defineSet("obs", "x"); - - _pdf = _ws->pdf("SumGausExpo"); - if (printLevel > 1) { - _ws->Print(); - } - generateData(nEvts); + std::unique_ptr data{generateData(gauss, {x, sigma})}; + randomizeParameters({mean}, 1337); + runFitBenchmark(state, gauss, *data); +} - if (printLevel > 1) { - _data->Print(); - } - } - RooAbsPdf &pdf() { return *_pdf; } - RooAbsData &data() { return *_data; } +/// Sum of a Gaussian and an exponential. +static void benchFitAddPdf(benchmark::State &state) +{ + auto changeMsgLevel = makeChangeMsgLvl(); + + RooRealVar x("x", "x", -1.5, 40.5); + RooRealVar mean("mean", "mean of gaussian", 20., -10, 30); + RooRealVar sigma("sigma", "width of gaussian", 4., 0.5, 10); + RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma); + RooRealVar c1("c1", "Decay constant", -0.5, -10, -0.001); + RooExponential ex{"Expo", "Exponential PDF", x, c1}; + RooRealVar fractionGaus("fractionGaus", "Fraction of Gauss component", 0.5, 0., 1.); + RooAddPdf addPdf("SumGausExpo", "Sum of Gaus and Exponential", {gauss, ex}, fractionGaus); + // to avoid a warning when computing the unnormalized RooAddPdf values + addPdf.fixCoefNormalization(x); + + std::unique_ptr data{generateData(addPdf, {x})}; + randomizeParameters({mean, sigma, c1, fractionGaus}, 1337); + runFitBenchmark(state, addPdf, *data); +} - void randomiseParameters(unsigned long seed = 0) - { - std::unique_ptr parameters{_pdf->getParameters(_data.get())}; - auto random = RooRandom::randomGenerator(); - if (seed != 0) - random->SetSeed(seed); - - for (auto param : *parameters) { - auto par = static_cast(param); - const double uni = random->Uniform(); - const double min = par->getMin(); - const double max = par->getMax(); - par->setVal(min + uni * (max - min)); - } - } +/// Product of three one-dimensional mixtures. +static void benchFitProdPdf(benchmark::State &state) +{ + auto changeMsgLevel = makeChangeMsgLvl(); - void generateData(size_t nEvts) - { - _data = std::unique_ptr(_pdf->generate(*_ws->set("obs"), nEvts)); - // allocate here output vector - _results = std::vector(nEvts); - } -}; + RooWorkspace w; + w.factory("Gaussian::g1(x[0,20],mx1[10,0,20],sx1[3,0.1,10])"); + w.factory("Gaussian::g2(x,mx2[10,0,20],sx2[1,0.1,10])"); + w.factory("Gaussian::g3(y[0,20],my[10,0,20],sy[2,0.1,10])"); + w.factory("Polynomial::pol(y,a[-0.01,-0.05,0.1])"); + w.factory("Gamma::gpdf(z[0,10],g[2,0.1,10],b[2,0.1,10],m0[-1,-20,0])"); + w.factory("SUM::pdfx(fx[0.2,0,1]*g1,g2)"); + w.factory("SUM::pdfy(fy[0.4,0.,1.]*g3,pol)"); + w.factory("PROD::model(pdfx,pdfy,gpdf)"); -static void benchFit(benchmark::State &state) -{ - TestModelPdf model(nEvents); - // We need to randomize the parameters so we don't start already at the - // minimum (the parameter values for which the dataset was created). - model.randomiseParameters(1337); - runFitBenchmark(state, model.pdf(), model.data()); + auto &pdf = *w.pdf("model"); + std::unique_ptr data{generateData(pdf, {*w.var("x"), *w.var("y"), *w.var("z")})}; + RooArgSet parameters; + pdf.getParameters(data->get(), parameters); + shiftParameters(parameters, 0.05); + runFitBenchmark(state, pdf, *data); } -static void benchModel(benchmark::State &state) +/// Extended sum of five components. +/// +/// Watch out with the result from this benchmark: if there are evaluation +/// errors during the minimization, the time differences will mostly come from +/// how the backends are handling those. This might or might not be wanted. +static void benchFitModel(benchmark::State &state) { - RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING); - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range + auto changeMsgLevel = makeChangeMsgLvl(); RooWorkspace w; - // w.factory("Voigtian::vpdf(x[-10,10],m0[0,-10,10],width[0.5,0,10],s0[0.5,0,10])"); - // w.factory("Gaussian::vpdf(x[0,20],m0[10,0,20],s0[0.5,0,10])"); w.factory("Gamma::vpdf(x[0,20],g[20,0.1,40],b[0.5,0.01,10],m0[-1,-20,0])"); w.factory("Gaussian::gpdf(x,m1[10,0,20],s1[2,0.1,10])"); w.factory("Gaussian::g2(x,m2[5,0,20],s2[0.3,0.01,10])"); w.factory("Gaussian::g3(x,m3[15,0,20],s3[0.4,0.01,10])"); w.factory("Polynomial::pol(x,a[-0.01,-0.05,0.1])"); - w.factory("SUM::model(ns1[0,1000000000]*vpdf,ns2[0,10000000000]*gpdf,ng2[0,1000000000]*g2,ng3[0,10000000000]*g3, " "npol[0,10000000000]*pol)"); - // set coeffients w.var("ns1")->setVal(0.2 * nEvents); w.var("ns2")->setVal(0.3 * nEvents); w.var("ng2")->setVal(0.1 * nEvents); w.var("ng3")->setVal(0.1 * nEvents); w.var("npol")->setVal(0.3 * nEvents); - auto &pdf = *(w.pdf("model")); - auto &x = *(w.var("x")); - - std::unique_ptr data{pdf.generate({x}, nEvents)}; - + auto &pdf = *w.pdf("model"); + std::unique_ptr data{generateData(pdf, {*w.var("x")})}; RooArgSet parameters; pdf.getParameters(data->get(), parameters); - // We need to change the parameters so we don't start already at the minimum - // (the parameter values for which the dataset was created). This time, we - // don't randomize them because the fit has many parameter and is not stable - // from many initial points. That's why they are only shifted away a little - // bit from their true value. - for (auto *param : static_range_cast(parameters)) { - double val = param->getVal(); - val = val + 0.01 * (param->getMax() - val); - param->setVal(val); - } - + shiftParameters(parameters, 0.01); runFitBenchmark(state, pdf, *data); } -auto const unit = benchmark::kMillisecond; +/// B decay with mixing and a delta function resolution model. +static void benchFitBDecayWithMixing(benchmark::State &state) +{ + auto changeMsgLevel = makeChangeMsgLvl(); + + BDecayVars v; + RooTruthModel tm("tm", "truth model", v.dt); + RooBMixDecay bmix("bmix", "decay", v.dt, v.mixState, v.tagFlav, v.tau, v.dm, v.w, v.dw, tm, + RooBMixDecay::DoubleSided); + v.generateAndRun(state, bmix); +} + +/// B decay convolved with a Gaussian resolution model. +static void benchFitBDecayGaussResolution(benchmark::State &state) +{ + auto changeMsgLevel = makeChangeMsgLvl(); + + BDecayVars v; + RooRealVar bias1("bias1", "bias1", 0, -1., 1.); + RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9); + RooGaussModel gm1("gm1", "gauss model 1", v.dt, bias1, sigma1); + RooDecay decay_gm1("decay_gm1", "decay", v.dt, v.tau, gm1, RooDecay::DoubleSided); + v.generateAndRun(state, decay_gm1); +} + +/// B decay convolved with a sum of two Gaussian resolution models. +static void benchFitBDecayDoubleGauss(benchmark::State &state) +{ + auto changeMsgLevel = makeChangeMsgLvl(); + + BDecayVars v; + RooRealVar bias1("bias1", "bias1", 0, -1., 1.); + RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9); + RooGaussModel gm1("gm1", "gauss model 1", v.dt, bias1, sigma1); + RooRealVar bias2("bias2", "bias2", 0, -1., 1.); + RooRealVar sigma2("sigma2", "sigma2", 5, 0.5, 10); + RooGaussModel gm2("gm2", "gauss model 2", v.dt, bias2, sigma2); + RooRealVar gm1frac("gm1frac", "fraction of gm1", 0.5, 0.0, 1.0); + RooAddModel gmsum("gmsum", "sum of gm1 and gm2", {gm1, gm2}, gm1frac); + RooDecay decay_gmsum("decay_gmsum", "decay", v.dt, v.tau, gmsum, RooDecay::DoubleSided); + v.generateAndRun(state, decay_gmsum); +} + +namespace { + +constexpr auto unit = benchmark::kMillisecond; + +constexpr auto Cpu = static_cast(RooFit::EvalBackend::Value::Cpu); +constexpr auto Cuda = static_cast(RooFit::EvalBackend::Value::Cuda); +constexpr auto Codegen = static_cast(RooFit::EvalBackend::Value::Codegen); + +} // namespace + +// The benchmark arguments are the evaluation backend and the number of +// threads, where zero means all threads of the implicit multi-threading pool. +#define BENCH_ENTRY(func, name, backend, nThreads) \ + BENCHMARK(func)->Unit(unit)->UseRealTime()->Name(name)->Args({backend, nThreads}); #ifdef DO_BENCH_ROOFIT_CUDA -// Only one iteration for CUDA benchmarks. Otherwise, the benchmark suit will -// repeat them many times are they appear to be very fast, consuming almost no -// CPU time. -#define CUDA_ONLY(code) code->Iterations(1) +#define BENCH_CUDA(func, model) BENCH_ENTRY(func, model "_FitCUDA", Cuda, 1) #else -#define CUDA_ONLY(code) +#define BENCH_CUDA(func, model) #endif -BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCPU")->Args({fitCpu}); -BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegen")->Args({fitCodegen}); -BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegenNoGrad")->Args({fitCodegenNoGrad}); -CUDA_ONLY(BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCUDA")->Args({fitCuda})); - -// No codegen variants for this model: normalizing the Gaussian over both x -// and sigma requires a multi-dimensional numeric integral, which the codegen -// backend doesn't support. -BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCPU")->Args({fitCpu}); -CUDA_ONLY(BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCUDA")->Args({fitCuda})); - -BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCPU")->Args({fitCpu}); -BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegen")->Args({fitCodegen}); -BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegenNoGrad")->Args({fitCodegenNoGrad}); -CUDA_ONLY(BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCUDA")->Args({fitCuda})); - -BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu}); -BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen}); -BENCHMARK(benchProdPdf) - ->Unit(unit) - ->Name("ProdPdf_FitCodegenNoGrad") - ->Unit(benchmark::kMillisecond) - ->Args({fitCodegenNoGrad}); -CUDA_ONLY(BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCUDA")->Unit(benchmark::kMillisecond)->Args({fitCuda})); - -// Watch out with the result from these benchmarks: if there are evaluation -// errors during the minimization, the time differences will mostly come from -// how the backends are handling those. This might or might not be wanted. -BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu}); -BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen}); -BENCHMARK(benchModel) - ->Unit(unit) - ->Name("FitModel_FitCodegenNoGrad") - ->Unit(benchmark::kMillisecond) - ->Args({fitCodegenNoGrad}); -CUDA_ONLY(BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCUDA")->Unit(benchmark::kMillisecond)->Args({fitCuda})); +// The CPU and CUDA backends support all models. +#define BENCH_BACKENDS(func, model) \ + BENCH_ENTRY(func, model "_FitCPU", Cpu, 1) \ + BENCH_ENTRY(func, model "_FitCPUMT", Cpu, 0) \ + BENCH_CUDA(func, model) + +// The codegen variant is only registered for models that the codegen backend +// supports, and only if ROOT was built with clad, which the gradient needs. +#ifdef R__HAS_CLAD +#define BENCH_BACKENDS_WITH_CODEGEN(func, model) \ + BENCH_BACKENDS(func, model) \ + BENCH_ENTRY(func, model "_FitCodegen", Codegen, 1) +#else +#define BENCH_BACKENDS_WITH_CODEGEN(func, model) BENCH_BACKENDS(func, model) +#endif -int main(int argc, char **argv) +BENCH_BACKENDS_WITH_CODEGEN(benchFitGauss, "Gaus") +// Normalizing the Gaussian over both x and sigma requires a multi-dimensional +// numeric integral, which the codegen backend doesn't support. +BENCH_BACKENDS(benchFitGaussXSigma, "GausXS") +BENCH_BACKENDS_WITH_CODEGEN(benchFitAddPdf, "AddPdf") +BENCH_BACKENDS_WITH_CODEGEN(benchFitProdPdf, "ProdPdf") +BENCH_BACKENDS_WITH_CODEGEN(benchFitModel, "FitModel") +// The codegen backend has no analytical integrals for RooBMixDecay, +// RooGaussModel and RooAddModel yet. +BENCH_BACKENDS(benchFitBDecayWithMixing, "BDecayWithMixing") +BENCH_BACKENDS(benchFitBDecayGaussResolution, "BDecayGaussResolution") +BENCH_BACKENDS(benchFitBDecayDoubleGauss, "BDecayDoubleGauss") + +namespace { + +/// The CPU model name from /proc/cpuinfo, or "unknown" elsewhere than Linux. +std::string cpuModelName() { - // Parse command line arguments - for (Int_t i = 1; i < argc; i++) { - std::string arg = argv[i]; - if (arg == "-v") { - std::cout << "---running in verbose mode" << std::endl; - printLevel = 1; - } else if (arg == "-vv") { - std::cout << "---running in very verbose mode" << std::endl; - printLevel = 2; - } - if (arg == "-n") { - nEvents = atoi(argv[++i]); + std::ifstream in("/proc/cpuinfo"); + std::string line; + while (std::getline(in, line)) { + if (line.rfind("model name", 0) == 0) { + auto pos = line.find(':'); + if (pos != std::string::npos) { + pos = line.find_first_not_of(" \t", pos + 1); + return pos == std::string::npos ? "unknown" : line.substr(pos); + } } } + return "unknown"; +} + +/// The name of the first NVIDIA GPU, as reported by nvidia-smi. +std::string gpuName() +{ + if (!gSystem->Which(gSystem->Getenv("PATH"), "nvidia-smi")) { + return "unknown"; + } + TString out = gSystem->GetFromPipe("nvidia-smi --query-gpu=name --format=csv,noheader 2> /dev/null"); + std::string name = out.Data(); + name = name.substr(0, name.find('\n')); + return name.empty() ? "unknown" : name; +} + +/// The RooBatchCompute library that RooFit loaded for the CPU backend, e.g. +/// "AVX2" or "GENERIC". It is chosen by the "RooFit.BatchCompute" setting of +/// the ROOT configuration (default: generic, "auto" for the fastest one that +/// the CPU supports) when the first likelihood is evaluated, so that is +/// triggered here with a minimal model. +std::string batchComputeLibrary() +{ + { + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + RooRealVar x("x", "x", -5, 5); + RooRealVar mean("mean", "mean", 0, -1, 1); + RooRealVar sigma("sigma", "sigma", 1, 0.5, 2); + RooGaussian gauss("gauss", "gauss", x, mean, sigma); + std::unique_ptr data{gauss.generate({x}, 10)}; + std::unique_ptr nll{gauss.createNLL(*data, RooFit::EvalBackend::Cpu())}; + nll->getVal(); + } + // The name is taken from the list of loaded libraries, because the + // RooBatchCompute interface that knows the architecture is not a public + // header. The CUDA library is not loaded at this point. + const std::string prefix = "libRooBatchCompute_"; + std::string libs = gSystem->GetLibraries(); + auto pos = libs.find(prefix); + if (pos == std::string::npos) { + return "unknown"; + } + pos += prefix.size(); + return libs.substr(pos, libs.find(".so", pos) - pos); +} - std::cout << "Using nEvents = " << nEvents << std::endl; +} // namespace + +int main(int argc, char **argv) +{ + // Take the "--events=N" option out of the arguments before Google + // Benchmark parses the rest. + std::vector args; + for (int i = 0; i < argc; ++i) { + const std::string arg = argv[i]; + const std::string prefix = "--events="; + if (arg.rfind(prefix, 0) == 0) { + nEvents = std::strtoull(arg.c_str() + prefix.size(), nullptr, 10); + if (nEvents == 0) { + std::cerr << "invalid number of events in \"" << arg << "\"" << std::endl; + return 1; + } + } else { + args.push_back(argv[i]); + } + } + argc = args.size(); + args.push_back(nullptr); + argv = args.data(); benchmark::Initialize(&argc, argv); + if (benchmark::ReportUnrecognizedArguments(argc, argv)) { + return 1; + } + + // Enable the implicit multi-threading up front, so that the process state + // is the same for all benchmarks no matter in which order they run. + ROOT::EnableImplicitMT(); + + benchmark::AddCustomContext("root_version", gROOT->GetVersion()); + benchmark::AddCustomContext("batchcompute_library", batchComputeLibrary()); +#ifdef R__HAS_VDT + benchmark::AddCustomContext("vdt", "yes"); +#else + benchmark::AddCustomContext("vdt", "no"); +#endif + benchmark::AddCustomContext("events", std::to_string(nEvents)); + benchmark::AddCustomContext("threads", std::to_string(ROOT::GetThreadPoolSize())); + benchmark::AddCustomContext("cpu", cpuModelName()); +#ifdef DO_BENCH_ROOFIT_CUDA + benchmark::AddCustomContext("gpu", gpuName()); +#endif + benchmark::RunSpecifiedBenchmarks(); + benchmark::Shutdown(); + return 0; } diff --git a/root/roofit/roofit/benchRooFitBackends_make_plot.py b/root/roofit/roofit/benchRooFitBackends_make_plot.py index 65ab1f2ed..a45d002c2 100644 --- a/root/roofit/roofit/benchRooFitBackends_make_plot.py +++ b/root/roofit/roofit/benchRooFitBackends_make_plot.py @@ -1,115 +1,224 @@ -# To get the data to produce this plot, run benchmark with: -# ./benchRooFitBackends -n 1000000 --benchmark_format=json --benchmark_out=benchRooFitBackends.json - -import ROOT -import numpy as np +# Plots the speedup of the multi-threaded CPU and the CUDA evaluation backends +# relative to the single-threaded CPU backend for all models of +# benchRooFitBackends, together with the information about the run. +# +# To get the data, run the benchmark with the RooBatchCompute library for the +# CPU architecture: the default is the generic one without vectorization, and +# "auto" picks the fastest library that the CPU supports. The setting is read +# from the ROOT configuration, e.g. a .rootrc file in the working directory: +# +# echo "RooFit.BatchCompute: auto" > .rootrc +# ./benchRooFitBackends --benchmark_filter='Fit(CPU|CPUMT|CUDA)/' \ +# --benchmark_repetitions=3 \ +# --benchmark_format=json --benchmark_out=benchRooFitBackends.json +# python3 benchRooFitBackends_make_plot.py benchRooFitBackends.json +# +# Add for example --events=1000000 to fit 1M events instead of the default +# 100k. The slower fits run for a single iteration, so the repetitions give a +# spread estimate; the plot uses the median over the repetitions when the +# file has one. The plot is written next to the JSON file as PDF and PNG, +# and the timings and speedups are printed as a table. + +import argparse import json +import os import sys -from collections import defaultdict - -ROOT.gROOT.SetBatch(True) -json_name = sys.argv[1] -base_name = json_name.split(".")[0] +parser = argparse.ArgumentParser(description="Plot the speedups of the RooFit evaluation backends from benchRooFitBackends.") +parser.add_argument("json_file", help="JSON output of benchRooFitBackends") +parser.add_argument("--output", help="base name of the PDF and PNG files (default: like the JSON file)") +parser.add_argument("--xmax", type=float, default=30.0, help="upper end of the speedup axis (default: 30)") +# Parsed before importing ROOT, so that --help stays instant and PyROOT never +# gets a look at sys.argv. +args = parser.parse_args() -f = open(json_name, "rb") -json_object = json.load(f)["benchmarks"] -f.close() +import ROOT # noqa: E402 -timing_dict = defaultdict(dict) +ROOT.gROOT.SetBatch(True) -labels_in_json = [] +json_name = args.json_file +base_name = args.output or os.path.splitext(json_name)[0] -for bm in json_object: - name, backend = bm["name"].split("/")[0].split("_Fit") - timing_dict[name][backend] = bm["real_time"] - if len(labels_in_json) == 0 or labels_in_json[-1] != name: - labels_in_json.append(name) +with open(json_name) as f: + results = json.load(f) -# labels_in_json = [ -# "Gaus", -# "GausXS", -# "AddPdf", -# "ProdPdf", -# "FitModel", -# ] - -backends = ["CUDA"] -reference_backend = "CPU" +context = results["context"] +# Timings per model and backend, in the order of the benchmarks in the file. +# With --benchmark_repetitions, the file has one entry per repetition and the +# aggregates (mean, median, ...); then the median is used. timings = {} +for bm in results["benchmarks"]: + if bm.get("run_type", "iteration") == "aggregate" and bm.get("aggregate_name") != "median": + continue + model, backend = bm["name"].split("/")[0].split("_Fit") + entry = timings.setdefault(model, {}) + if bm.get("run_type") == "aggregate" or backend not in entry: + entry[backend] = bm["real_time"] + +# The thread count of the multi-threaded CPU backend: from the context, or +# from the counter of the benchmark for older files. +num_threads = context.get("threads") +if num_threads is None: + for bm in results["benchmarks"]: + if "_FitCPUMT" in bm["name"] and bm.get("num_threads"): + num_threads = int(bm["num_threads"]) + break +num_threads = int(num_threads) if num_threads else None -print(timing_dict) - -for b in backends + [reference_backend]: - l = list() - - for label in labels_in_json: - l.append(timing_dict[label][b]) - - timings[b] = np.array(l) - -# normalize to ratios -for b in backends: - timings[b] = timings[reference_backend] / timings[b] - -print(timings) +reference_backend = "CPU" +# Bars from top to bottom per model, with the legend entries. +backend_labels = {"CPUMT": "CPU, multi-threaded", "CUDA": "CUDA"} # CPUMT gets the thread count in the legend +backends = [b for b in backend_labels if any(b in t for t in timings.values())] + +models = [m for m in timings if reference_backend in timings[m]] +if not models: + sys.exit(f"no {reference_backend} results in {json_name}") +if not backends: + sys.exit(f"no {' or '.join(backend_labels)} results in {json_name}") + + +def speedup(model, backend): + """Speedup of a backend relative to the reference, or None if it was not run.""" + time = timings[model].get(backend) + return timings[model][reference_backend] / time if time else None + + +speedups = {b: [speedup(m, b) for m in models] for b in backends} + +# Print the timings and speedups as a table. +columns = [reference_backend] + backends +print(f"{'model':24s}" + "".join(f"{c + ' [ms]':>14s}" for c in columns) + "".join(f"{c + ' x':>10s}" for c in backends)) +for model in models: + line = f"{model:24s}" + for c in columns: + time = timings[model].get(c) + line += f"{time:14.1f}" if time else f"{'-':>14s}" + for b in backends: + value = speedup(model, b) + line += f"{value:10.1f}" if value else f"{'-':>10s}" + print(line) labels_dict = { "Gaus": "Gaussian", "GausXS": "#splitline{Gaussian}{(with #sigma as observable)}", "AddPdf": "#splitline{Gaussian}{plus exponential}", - "ProdPdf": "#splitline{Mixture 1:}{#splitline{(gauss + gauss)}{#splitline{x (gauss + poly)}{ x gamma}}}", - "FitModel": "#splitline{#splitline{Mixture 2:}{gamma + gamma}}{ + gauss + gauss + poly}", + "ProdPdf": "#splitline{(gauss + gauss)}{#splitline{#times (gauss + poly)}{#times gamma}}", + "FitModel": "#splitline{gamma + 3 gauss}{+ poly, extended}", "BDecayWithMixing": "#splitline{RooBMixDecay}{perfect resolution}", - "BDecayGaussResolution": "#splitline{RooBMixDecay}{Gaussian res.}", - "BDecayDoubleGauss": "#splitline{RooBMixDecay}{Double-Gaussian res.}", + "BDecayGaussResolution": "#splitline{RooDecay}{Gaussian res.}", + "BDecayDoubleGauss": "#splitline{RooDecay}{double-Gaussian res.}", } -labels = [labels_dict[l] for l in labels_in_json] - -nx = len(timings[reference_backend]) -nb = len(backends) - -# cb = ROOT.TCanvas("cb","cb",900,400) -cb = ROOT.TCanvas("cb", "cb", int(180 * len(labels) * 2), int(400 * 2)) -cb.SetGrid() -cb.SetBottomMargin(0.16) - -# ROOT.gROOT.GetColor(9).SetRGB(12.2 / 100, 46.7 / 100, 70.6 / 100) -# ROOT.gROOT.GetColor(46).SetRGB(100.0 / 100, 49.8 / 100, 5.5 / 100) -# ROOT.gROOT.GetColor(9).SetRGB(17.3 / 100, 62.7 / 100, 17.3 / 100) - -y_max = 1 -for b in backends: - y_max = max(y_max, np.max(timings[b]) + 3) +# Run information, shown as a two-column table above the plot. +events = context.get("events", "") +events = f"{int(events):,}" if events.isdigit() else "unknown" +info_rows = [ + ("CPU", context.get("cpu", "unknown")), + ("GPU", context.get("gpu", "none")), + ("Events per fit", events), + ("CPU threads", (str(num_threads) + " (multi-threaded CPU backend)") if num_threads else "unknown"), + ("Software", "ROOT " + context.get("root_version", "unknown") + ", RooBatchCompute " + + context.get("batchcompute_library", "unknown") + ", VDT: " + context.get("vdt", "unknown")), + ("Run", context.get("date", "")[:10] + " on " + context.get("host_name", "unknown")), +] + +n_models = len(models) +n_backends = len(backends) +# The range of the speedup axis is fixed, such that plots from different runs +# can be compared visually. +x_max = args.xmax + +# Horizontal bars: the models are on the vertical axis (the "x" axis of the +# histogram drawn with the "hbar" option) and the speedup on the horizontal +# one. The first model is drawn at the top. The top margin holds the title and +# the run information. +canvas = ROOT.TCanvas("canvas", "canvas", 1400, 420 + 120 * n_models) +left_margin = 0.25 +top_margin = 0.26 +canvas.SetLeftMargin(left_margin) +canvas.SetRightMargin(0.03) +canvas.SetTopMargin(top_margin) +canvas.SetBottomMargin(0.10) +canvas.SetGridx() ROOT.gStyle.SetHistMinimumZero() - -h1b = ROOT.TH1F( - "h1b", "RooFit: speedup in benchmark fits with the CUDA backend relative to the CPU backend (1 million events)", nx, 0, nx -) -h1b.SetFillColor(9) -h1b.SetBarWidth(0.8 / nb) -h1b.SetBarOffset(0.1) -h1b.SetStats(0) -h1b.SetMinimum(-0.1) -h1b.SetMaximum(y_max) - -h1b.GetYaxis().SetTitle("Relative Speedup") -h1b.GetYaxis().SetTitleSize(0.045) -h1b.GetXaxis().SetLabelSize(0.055) -# h1b.GetXaxis().SetLabelOffset(0.1) - -for i in range(nx): - h1b.SetBinContent(i + 1, timings["CUDA"][i]) - h1b.GetXaxis().SetBinLabel(i + 1, labels[i]) - -h1b.Draw("b") - -legend = ROOT.TLegend(0.1, 0.7, 0.46, 0.9) -legend.AddEntry(h1b, "CUDA relative to CPU (single thread)", "f") +colors = {"CPUMT": ROOT.kAzure + 1, "CUDA": ROOT.kGreen + 2} + +hists = [] +bar_width = 0.8 / n_backends +for i_backend, backend in enumerate(backends): + h = ROOT.TH1F(f"h_{backend}", "", n_models, 0, n_models) + for i_model, model in enumerate(models): + # Reverse the order such that the first model is at the top. + h.SetBinContent(n_models - i_model, speedups[backend][i_model] or 0.0) + h.GetXaxis().SetBinLabel(n_models - i_model, labels_dict.get(model, model)) + h.SetFillColor(colors.get(backend, ROOT.kGray)) + h.SetBarWidth(bar_width) + h.SetBarOffset(0.1 + (n_backends - 1 - i_backend) * bar_width) + h.SetStats(False) + h.SetMinimum(0.0) + h.SetMaximum(x_max) + hists.append(h) + +first = hists[0] +first.GetXaxis().SetLabelSize(0.038) +first.GetYaxis().SetTitle("Migrad speedup relative to the single-threaded CPU backend") +first.GetYaxis().SetLabelSize(0.038) +first.GetYaxis().SetTitleSize(0.038) +first.GetYaxis().SetTitleOffset(1.1) +first.Draw("hbar") +for h in hists[1:]: + h.Draw("hbar same") + +# Print the speedup next to each bar. Bars beyond the axis range are cut, and +# their value is printed at the end of the axis. +texts = [] +for i_backend, backend in enumerate(backends): + for i_model in range(n_models): + value = speedups[backend][i_model] + if value is None: + continue + y = n_models - i_model - 1 + 0.1 + (n_backends - 1 - i_backend + 0.5) * bar_width + text = ROOT.TLatex(min(value, 0.9 * x_max) + 0.01 * x_max, y, f"{value:.1f}#times") + text.SetTextAlign(12) + text.SetTextSize(0.034) + text.Draw() + texts.append(text) + +# The legend goes into the free space above the model labels, next to the +# run information. +legend = ROOT.TLegend(0.01, 0.78, 0.23, 0.90) +legend.SetBorderSize(0) +legend.SetMargin(0.16) +legend.SetFillStyle(0) +legend.SetTextSize(0.028) +for h, backend in zip(hists, backends): + label = backend_labels[backend] + if backend == "CPUMT" and num_threads: + label = f"CPU, {num_threads} threads" + legend.AddEntry(h, label, "f") legend.Draw() -cb.SaveAs(base_name + ".png") -cb.SaveAs(base_name + ".pdf") +title = ROOT.TLatex(0.5, 0.965, "RooFit evaluation backends in unbinned fits (Migrad, strategy 0)") +title.SetNDC() +title.SetTextAlign(22) +title.SetTextSize(0.036) +title.Draw() + +# The two columns of the run information, aligned with the frame. +info_texts = [] +row_height = 0.027 +for i_row, (key, value) in enumerate(info_rows): + y = 0.925 - (i_row + 1) * row_height + for x, text, font in ((left_margin, key, 62), (left_margin + 0.15, value, 42)): + latex = ROOT.TLatex(x, y, text) + latex.SetNDC() + latex.SetTextAlign(12) + latex.SetTextFont(font) + latex.SetTextSize(0.025) + latex.Draw() + info_texts.append(latex) + +canvas.SaveAs(base_name + ".pdf") +canvas.SaveAs(base_name + ".png")