Benchmark: time per toy

[mmikhasenko]

A benchmark from Morefit preprint,

Generate data + fit with Minuit,

Result of different frameworks:

Here is a first, naive test of Julia using Distributions.jl + Minuit2.jl

using Plots
using FHist
using Optim
using QuadGK
using Minuit2
using Parameters
using DataFrames
using Distributions
using ComponentArrays

function total_model(pars, a, b)
	@unpack log_fb = pars
	prior = [1, exp(log_fb)] ./ (1+exp(log_fb))
	MixtureModel(
		[signal_model(pars.sig, a, b), bgd_model(pars.bgd, a, b)],
		prior)
end

function signal_model(pars, a, b)
	@unpack μ, σ = pars
	return truncated(Normal(μ, σ), a, b)
end

function bgd_model(pars, a, b)
	@unpack τ = pars
	return truncated(Exponential(τ), a, b)
end

const gen_pars = ComponentArray(
	sig=(μ=5.28, σ=0.06),
	bgd=(τ=1.0,),
	log_fb = 1/(1/(1-0.7)-1)
)

const x_range = (5.0, 7.0)
model = total_model(gen_pars, x_range...)

# plot(x->pdf(model,x), x_range..., )

With Minuit2

function fit_nll_minuit(build_model, data; init_pars)
	ax = getaxes(init_pars)
	function objective(x)
		pars = (;
			sig=(μ=x[1], σ=x[2]),
			bgd=(τ=x[3],),
			log_fb = x[4])
		nll(build_model(pars), data)
	end
	m = Minuit(objective, init_pars; tolerance=1e-6)
	migrad!(m)
	#
	best_pars = ComponentArray(collect(m.values), getaxes(init_pars))
	best_model = build_model(best_pars)
	return (; best_pars, best_model, fit_res = m)
end

@btime a_toy_minuit(1_000).fit_res  #   18.706 ms (896360 allocations: 17.87 MiB)
@btime a_toy_minuit(10_000).fit_res  #   126.391 ms (9386461 allocations: 186.56 MiB)
@btime a_toy_minuit(100_000).fit_res  #   1.505 s (105806691 allocations: 2.04 GiB)

# roughly 100 calls, 7 iterations

With Optim, Nelder-Mead

function fit_nll_optim(build_model, data; init_pars)
	objective(pars) = nll(build_model(pars), data)
	fit_res = optimize(objective, init_pars)
	best_pars = fit_res.minimizer
	best_model = build_model(best_pars)
	return (; best_pars, best_model, fit_res)
end

function a_toy_optim(n)
	data = rand(model, n)

	# fit
	build_model(pars) = total_model(pars, x_range...)
	(; best_pars, best_model, fit_res) =
		fit_nll_optim(build_model, data; init_pars = gen_pars)
	
	(; fit_res.iterations, fit_res.f_calls, fit_res.g_calls, fit_res.minimum)
end

@btime a_toy_optim(1000)  # 31.942 ms (2442839 allocations: 48.26 MiB)
@btime a_toy_optim(10_000) # 362.661 ms (28705762 allocations: 570.04 MiB)
@btime a_toy_optim(100_000) # 3.806 s (294205916 allocations: 5.67 GiB)

# roughly 100 calls, 

System

# system
Darwin Kernel Version 24.5.0: Tue Apr 22 19:54:33 PDT 2025; root:xnu-11417.121.6~2/RELEASE_ARM64_T8122

# julia
  [6e4b80f9] BenchmarkTools v1.6.0
  [b0b7db55] ComponentArrays v0.15.27
  [a93c6f00] DataFrames v1.7.0
  [31c24e10] Distributions v0.25.120
  [68837c9b] FHist v0.11.11
⌅ [37821647] Minuit2 v0.2.1
⌅ [429524aa] Optim v1.12.0
  [d96e819e] Parameters v0.12.3
  [91a5bcdd] Plots v1.40.13
  [1fd47b50] QuadGK v2.11.2
  [44cfe95a] Pkg v1.11.0

[peremato]

@mmikhasenko I plan to work on this benchmark asap