diff --git a/root/roofit/roofit/benchRooFitBackends.cxx b/root/roofit/roofit/benchRooFitBackends.cxx index eaa90481..6b633e0e 100644 --- a/root/roofit/roofit/benchRooFitBackends.cxx +++ b/root/roofit/roofit/benchRooFitBackends.cxx @@ -14,7 +14,7 @@ // 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 +// created outside of the timed loop and is offset by its initial value. 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 @@ -139,7 +139,12 @@ void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data) RooArgSet paramsInitial; params.snapshot(paramsInitial); - std::unique_ptr nll{pdf.createNLL(data, EvalBackend(evalBackend), NumCPU(nThreads))}; + // The likelihood is offset by its initial value, like it should be done for + // large datasets: without it, the rounding errors of the large absolute NLL + // values make Migrad fail to reach the required EDM from a few million + // events on. + std::unique_ptr nll{ + pdf.createNLL(data, EvalBackend(evalBackend), NumCPU(nThreads), Offset("initial"))}; // Every fit gets a fresh minimizer: the evaluation error recovery state of // the minimizer function (maximum FCN value, number of invalid @@ -285,27 +290,30 @@ static void benchFitProdPdf(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. +/// The Gamma component is a shoulder at low x, distinct from the broad +/// Gaussian, and the yield ranges scale with the number of events. Otherwise +/// the yields of two nearly identical components are degenerate and Migrad +/// often fails to converge, in particular with many events. static void benchFitModel(benchmark::State &state) { auto changeMsgLevel = makeChangeMsgLvl(); RooWorkspace w; - w.factory("Gamma::vpdf(x[0,20],g[20,0.1,40],b[0.5,0.01,10],m0[-1,-20,0])"); + w.factory("Gamma::vpdf(x[0,20],g[2,0.1,40],b[2,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)"); - - 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); + w.factory("SUM::model(ns1[0,1]*vpdf,ns2[0,1]*gpdf,ng2[0,1]*g2,ng3[0,1]*g3,npol[0,1]*pol)"); + + // The yields are set to the expected numbers of events, with ranges up to + // ten times the total number of events. + const std::pair yields[] = { + {"ns1", 0.2}, {"ns2", 0.3}, {"ng2", 0.1}, {"ng3", 0.1}, {"npol", 0.3}}; + for (auto const &[name, fraction] : yields) { + w.var(name)->setRange(0, 10. * nEvents); + w.var(name)->setVal(fraction * nEvents); + } auto &pdf = *w.pdf("model"); std::unique_ptr data{generateData(pdf, {*w.var("x")})}; diff --git a/root/roofit/roofit/benchRooFitBackends_make_plot.py b/root/roofit/roofit/benchRooFitBackends_make_plot.py index a45d002c..2ad73ba9 100644 --- a/root/roofit/roofit/benchRooFitBackends_make_plot.py +++ b/root/roofit/roofit/benchRooFitBackends_make_plot.py @@ -1,6 +1,10 @@ -# 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. +# Plots the cost of the RooFit evaluation backends for all models of +# benchRooFitBackends: the median wall time of a fit divided by the number of +# likelihood evaluations and the number of events, so that models with very +# different fit times and runs with different numbers of events can be shown +# on one axis. The speedups between the backends are distances on the +# logarithmic axis. A second run can be overlaid with --compare, e.g. to +# compare machines; its markers are hollow. # # To get the data, run the benchmark with the RooBatchCompute library for the # CPU architecture: the default is the generic one without vectorization, and @@ -11,23 +15,30 @@ # ./benchRooFitBackends --benchmark_filter='Fit(CPU|CPUMT|CUDA)/' \ # --benchmark_repetitions=3 \ # --benchmark_format=json --benchmark_out=benchRooFitBackends.json -# python3 benchRooFitBackends_make_plot.py benchRooFitBackends.json +# python3 benchRooFitBackends_make_plot.py benchRooFitBackends.json \ +# [--compare other_machine.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. +# and the timings are printed as a table. import argparse import json import os +import re import sys -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 = argparse.ArgumentParser( + description="Plot the cost per likelihood evaluation and event of the RooFit evaluation backends " + "from benchRooFitBackends." +) +parser.add_argument("json_file", help="JSON output of benchRooFitBackends (filled markers)") +parser.add_argument("--compare", metavar="JSON_FILE", help="second run to overlay (hollow markers)") 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)") +parser.add_argument("--xmin", type=float, default=0.1, help="lower end of the axis in ns (default: 0.1)") +parser.add_argument("--xmax", type=float, default=1000.0, help="upper end of the axis in ns (default: 1000)") # Parsed before importing ROOT, so that --help stays instant and PyROOT never # gets a look at sys.argv. args = parser.parse_args() @@ -36,68 +47,84 @@ ROOT.gROOT.SetBatch(True) -json_name = args.json_file -base_name = args.output or os.path.splitext(json_name)[0] - -with open(json_name) as f: - results = json.load(f) - -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: +base_name = args.output or os.path.splitext(args.json_file)[0] + +# Backends in the order of the legend, with marker styles (filled, hollow) and colors. +backend_styles = { + "CPU": ("CPU, 1 thread", 21, 25, ROOT.kGray + 2), + "CPUMT": ("CPU, all threads", 22, 26, ROOT.kAzure + 1), + "CUDA": ("CUDA", 20, 24, ROOT.kGreen + 2), +} + + +def short_cpu_name(name): + """E.g. "AMD Ryzen 9 5950X 16-Core Processor" -> "Ryzen 9 5950X".""" + name = re.sub(r"\((R|TM)\)", "", name) + name = re.sub(r"\s+\d+-Core Processor|\s+CPU\s+@.*|^(AMD|Intel)\s+", "", name) + return name.strip() + + +def load_run(path): + """Read one JSON file. For every model and backend, the fit time is the + median over the repetitions if the file has one, otherwise the first + iteration entry. Returns the context and {model: {backend: entry}}.""" + with open(path) as f: + results = json.load(f) + context = results["context"] + events = int(context["events"]) if context.get("events", "").isdigit() else None + if events is None: + sys.exit(f"{path}: the number of events is missing in the benchmark context") + data = {} 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 - -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) + is_aggregate = bm.get("run_type", "iteration") == "aggregate" + if is_aggregate and bm.get("aggregate_name") != "median": + continue + model, backend = bm["name"].split("/")[0].split("_Fit") + if backend not in backend_styles: + continue + entry = data.setdefault(model, {}) + if backend in entry and not is_aggregate: + continue + evals = bm.get("nll_evals") + if not evals: + sys.exit(f"{path}: no nll_evals counter for {bm['name']}; the file is from an older benchmark version") + entry[backend] = { + "time_ms": bm["real_time"], + "evals": evals, + "converged": bm.get("migrad_status", 0) == 0, + # Fit time per likelihood evaluation and event in ns. + "ns": bm["real_time"] * 1e6 / (evals * events), + } + return {"context": context, "events": events, "data": data, "path": path} + + +runs = [load_run(args.json_file)] +if args.compare: + runs.append(load_run(args.compare)) +compare = len(runs) > 1 + +models = list(runs[0]["data"]) +backends = [b for b in backend_styles if any(b in r["data"].get(m, {}) for r in runs for m in models)] +if not models or not backends: + sys.exit(f"no results in {args.json_file}") + +# Print the numbers as a table, one block per run. +for run in runs: + print(f"{run['path']}: {run['context'].get('cpu', 'unknown')}, {run['events']:,} events") + print(f"{'model':24s}" + "".join(f"{b + ' [ms]':>14s}{'evals':>7s}{'ns':>9s}" for b in backends)) + for model in models: + line = f"{model:24s}" + for b in backends: + e = run["data"].get(model, {}).get(b) + if e is None: + line += f"{'-':>14s}{'-':>7s}{'-':>9s}" + else: + flag = "" if e["converged"] else "*" + line += f"{e['time_ms']:14.1f}{int(e['evals']):7d}{e['ns']:8.3g}{flag:1s}" + print(line) + print() +if any(not e["converged"] for r in runs for d in r["data"].values() for e in d.values()): + print("* Migrad did not converge (status != 0), so the backends may have done different work.") labels_dict = { "Gaus": "Gaussian", @@ -110,115 +137,200 @@ def speedup(model, backend): "BDecayDoubleGauss": "#splitline{RooDecay}{double-Gaussian res.}", } -# 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")), -] + +def info_values(run): + """Values of the run information table for one run. The values are + shortened when two runs are shown next to each other.""" + c = run["context"] + cpu = c.get("cpu", "unknown") + gpu = c.get("gpu", "none") + software = f"{c.get('root_version', '?')}, {c.get('batchcompute_library', '?')}, VDT: {c.get('vdt', '?')}" + run_info = c.get("date", "")[:10] + (", " if compare else " on ") + c.get("host_name", "unknown") + if compare: + cpu = short_cpu_name(cpu) + gpu = re.sub(r"^NVIDIA\s+", "", gpu) + else: + software = "ROOT " + software.replace(", ", ", RooBatchCompute ", 1) + threads = c.get("threads", "unknown") + return [cpu, gpu, f"{run['events']:,}", threads if compare else f"{threads} (multi-threaded CPU backend)", + software, run_info] + + +info_labels = ["CPU", "GPU", "Events per fit", "CPU threads", "ROOT" if compare else "Software", "Run"] 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) +x_min, x_max = args.xmin, args.xmax + +# The layout is defined in pixels and converted to the fractions of the canvas +# height that ROOT uses, such that the text has the same size no matter how +# many rows there are. With two runs, the rows are higher to make room for two +# rows of markers. +top_px = 386 if compare else 359 +bottom_px = 138 +row_px = 180 if compare else 120 +height = top_px + bottom_px + row_px * n_models + + +def size(px): + """Text size in pixels as a fraction of the canvas height.""" + return px / height + + +def y_ndc(px_from_top): + return 1.0 - px_from_top / height + + +canvas = ROOT.TCanvas("canvas", "canvas", 1400, height) 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.SetRightMargin(0.045) +canvas.SetTopMargin(top_px / height) +canvas.SetBottomMargin(bottom_px / height) +canvas.SetLogx() canvas.SetGridx() -ROOT.gStyle.SetHistMinimumZero() -colors = {"CPUMT": ROOT.kAzure + 1, "CUDA": ROOT.kGreen + 2} +# Empty frame with one row per model, the first model at the top. The axis +# range is fixed, such that plots from different runs can be compared visually. +frame = ROOT.TH2F("frame", "", 1, x_min, x_max, n_models, 0, n_models) +frame.SetStats(False) +for i_model, model in enumerate(models): + frame.GetYaxis().SetBinLabel(n_models - i_model, labels_dict.get(model, model)) +frame.GetYaxis().SetLabelSize(size(52)) +frame.GetYaxis().SetTickLength(0) +frame.GetXaxis().SetTitle("Fit time / (likelihood evaluations #times events) [ns]") +frame.GetXaxis().SetLabelSize(size(52)) +frame.GetXaxis().SetTitleSize(size(52)) +frame.GetXaxis().SetTitleOffset(1.1) +frame.GetXaxis().SetNoExponent(True) +frame.Draw("AXIS") + +# Thin separators between the models. +keep = [] +for i in range(1, n_models): + line = ROOT.TLine() + line.SetLineColor(ROOT.kGray) + line.SetLineStyle(3) + line.DrawLine(x_min, i, x_max, i) + keep.append(line) + + +def value_text(value): + return f"{value:.2g}" if value < 100 else f"{value:.0f}" -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) + +# Vertical offsets of the marker rows and of the value labels within a model +# row, in units of rows. +marker_offsets = [26 / row_px, -26 / row_px] if compare else [0.0] +text_offset = 14 / row_px +value_size = size(28 if compare else 33) +for i_run, run in enumerate(runs): 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: + entries = run["data"].get(model, {}) + y = n_models - i_model - 0.5 + marker_offsets[i_run] + values = [entries[b]["ns"] for b in backends if b in entries] + if not values: 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) + # A line spanning the backends of this model, so the gaps between them + # (the speedups) are visible as lengths. + span = ROOT.TGraph(2) + span.SetPoint(0, min(values), y) + span.SetPoint(1, max(values), y) + span.SetLineColor(ROOT.kGray + 1) + span.SetLineWidth(2) + span.Draw("L") + keep.append(span) + for b in backends: + e = entries.get(b) + if e is None: + continue + _, filled, hollow, color = backend_styles[b] + marker = ROOT.TGraph(1) + marker.SetPoint(0, e["ns"], y) + marker.SetMarkerStyle(filled if i_run == 0 else hollow) + marker.SetMarkerColor(color) + marker.SetMarkerSize(2.6) + marker.Draw("P") + keep.append(marker) + # The value above the markers of the first run and below the + # markers of the second one. + above = i_run == 0 + text = ROOT.TLatex(e["ns"], y + (text_offset if above else -text_offset), + value_text(e["ns"]) + ("" if e["converged"] else "*")) + text.SetTextAlign(21 if above else 23) + text.SetTextSize(value_size) + text.SetTextColor(color if b != "CPU" else ROOT.kBlack) + text.Draw() + keep.append(text) + +# The legend for the backends goes into the free space above the model labels. +legend_top = y_ndc(110) +legend = ROOT.TLegend(0.01, legend_top - size(55) * len(backends), 0.23, legend_top) 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.SetTextSize(size(36)) +legend_markers = [] +for b in backends: + label, filled, _, color = backend_styles[b] + sample = ROOT.TGraph(1) + sample.SetMarkerStyle(filled) + sample.SetMarkerColor(color) + sample.SetMarkerSize(2.2) + legend.AddEntry(sample, label, "p") + legend_markers.append(sample) legend.Draw() -title = ROOT.TLatex(0.5, 0.965, "RooFit evaluation backends in unbinned fits (Migrad, strategy 0)") +if any(not e["converged"] for r in runs for d in r["data"].values() for e in d.values()): + note = ROOT.TLatex(0.02, legend_top - size(55) * len(backends) - size(25), "* Migrad did not converge") + note.SetNDC() + note.SetTextAlign(12) + note.SetTextSize(size(26)) + note.Draw() + keep.append(note) + +title = ROOT.TLatex(0.5, y_ndc(48), "RooFit evaluation backends in unbinned fits (Migrad, strategy 0)") title.SetNDC() title.SetTextAlign(22) -title.SetTextSize(0.036) +title.SetTextSize(size(50)) 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)): +# The run information as a table above the frame: one column of labels and +# one column of values per run. With two runs, the column headers show the +# marker style of each run. +info_size = size(30 if compare else 34.5) +row_height = size(34.5 if compare else 39) +label_width = 0.13 if compare else 0.15 +column_width = 0.29 +y_top = y_ndc(103) +first_row = 0 +if compare: + first_row = 1 + for i_run, run in enumerate(runs): + x = left_margin + label_width + i_run * column_width + sample = ROOT.TMarker(x + 0.008, y_top - row_height, 20 if i_run == 0 else 24) + sample.SetNDC() + sample.SetMarkerSize(1.8) + sample.Draw() + keep.append(sample) + header = ROOT.TLatex(x + 0.022, y_top - row_height, "filled markers" if i_run == 0 else "hollow markers") + header.SetNDC() + header.SetTextAlign(12) + header.SetTextFont(62) + header.SetTextSize(info_size) + header.Draw() + keep.append(header) +columns = [info_values(run) for run in runs] +for i_row, label in enumerate(info_labels): + y = y_top - (first_row + i_row + 1) * row_height + cells = [(left_margin, label, 62)] + cells += [(left_margin + label_width + i_run * column_width, columns[i_run][i_row], 42) for i_run in range(len(runs))] + for x, text, font in cells: latex = ROOT.TLatex(x, y, text) latex.SetNDC() latex.SetTextAlign(12) latex.SetTextFont(font) - latex.SetTextSize(0.025) + latex.SetTextSize(info_size) latex.Draw() - info_texts.append(latex) + keep.append(latex) canvas.SaveAs(base_name + ".pdf") canvas.SaveAs(base_name + ".png")