[FR] `PredictNumItersNeeded()` 1.4 correction factor

[ruben-laso]

Problem description
In the function PredictNumItersNeeded() there is this 1.4 correction factor.
This causes the time running the experiment to exceed by ~40% the time specified by --benchmark_min_time.
Of course, --benchmark_min_time denotes the minimum amount of time to run the benchmark, but an overrun of 40% seems excessive.
This is particularly relevant in supercomputers, where CPU time is expensive.

• Why is this estimation done, instead of stopping the iterations when the accumulated "iteration time" exceeds the target time?
• Is there a reason for selecting 1.4 as a correction factor?
• In cases where the execution times are not stable, could this prediction be wrong by a large margin?

Suggested solution
I suggest either removing the correction factor or making it configurable (with a default value of 1.0).

Example
As shown in the following output (executed with --benchmark_min_time=1s) the real execution time is ~1.4s: $7039 \times 198481 = 1397107759$, $8775 \times 160984 = 1412634600$, ...

------------------------------------------------------------------------------------------------------------------------
Benchmark                                                              Time             CPU   Iterations UserCounters...
------------------------------------------------------------------------------------------------------------------------
GNU-TBB/std::adjacent_difference/double/1024/manual_time            7039 ns         7022 ns       198481 bytes_per_second=6.50968Gi/s
GNU-TBB/std::adjacent_find/double/1024/manual_time                  8775 ns         8702 ns       160984 bytes_per_second=2.61075Gi/s
GNU-TBB/std::all_of/double/1024/manual_time                         7704 ns         7529 ns       199585 bytes_per_second=2.97387Gi/s
GNU-TBB/std::any_of/double/1024/manual_time                         4707 ns         4625 ns       301878 bytes_per_second=4.86754Gi/s

When executing the same code with --benchmark_min_time=0.71s ($1/1.4 \simeq 0.71$), the execution times are much closer to 1s: $7681 \times 134530 = 1033324930$, $9265 \times 103410 = 958093650$, ...

------------------------------------------------------------------------------------------------------------------------
Benchmark                                                              Time             CPU   Iterations UserCounters...
------------------------------------------------------------------------------------------------------------------------
GNU-TBB/std::adjacent_difference/double/1024/manual_time            7681 ns         7615 ns       134530 bytes_per_second=5.96565Gi/s
GNU-TBB/std::adjacent_find/double/1024/manual_time                  9265 ns         9175 ns       103410 bytes_per_second=2.47288Gi/s
GNU-TBB/std::all_of/double/1024/manual_time                         7996 ns         7572 ns       110333 bytes_per_second=2.86519Gi/s
GNU-TBB/std::any_of/double/1024/manual_time                         4656 ns         4680 ns       192865 bytes_per_second=4.92025Gi/s

[dmah42]

the multiplier is applied when we've run for less than min time to figure out how many iterations we need to get over the minimum time. we could go with a smaller multiplier, but it will take longer to get up to the minimum time, which (i believe) nets out at the same amount of computation at the end of the day.

[ruben-laso]

The confusing point is that the library will (internally) aim to run the benchmark for "min_time" + 40%.
Users might expect the benchmark to run just enough iterations to exceed the 'min_time' threshold.

What would happen if the estimation is wrong? Imagine an extreme situation where the first iteration(s) have very low performance due to a cold cache, but it is very fast after that.

What would be the implications of setting this number to 1.0? Would the function be called several times until the "min_time" is reached?

[dmah42]

if it was set to 1.0 then we would never hit mintime.

is it true we'd run for min_time + 40%? i think if the first run we try runs longer than mintime then we'd be done. the only time the multiplier would kick in is if an attempt runs for shorter than that time, at which point we'd increase the time of the next run by 40%. the only time we'd run for a time mintime * 1.4 is if the initial run is almost-but-not-quite mintime.

[LebedevRI]

As

benchmark/src/benchmark_runner.cc

Line 317 in 08fdf6e

double multiplier = GetMinTimeToApply() * 1.4 / std::max(i.seconds, 1e-9);

notes, we:

1. Run N (start=1) iteration
2. See if that is >= min time, if it is then we're done
3. If not, is it less-or-equal than 10% of min time? in which case we goto step 1 with N*=10
4. But if it is greater than 10% of min time, then we run somewhere between 10xN and 14xN iterations

So yes, we do need that fudge factor, and yes, we clearly can run for 1.4x of the min time.
In reality even more because of those initial runs with exponentially larger iteration times.

I'm guessing the factor of 1.4 is "derived" from statistical distribution of the iteration timings.
If we underestimate it, and the actual iteration times happen to be faster than the previous ones,
we may fall slightly short of the min time (say 0.9x of min time), we'll need to re-run again,
and that will end up being ~1.9x overshoot, which more wasteful than 1.4x overshoot.

[ruben-laso]

According to my calculations, in step 4 you would run somewhere between $1.4\times N$ and $14\times N$.

Since GetMinTimeToApply() $>$ i.seconds, we always have GetMinTimeToApply() / i.seconds $>1$.
Given that i.seconds / GetMinTimeToApply() $> 0.1$, we have GetMinTimeToApply() / i.seconds $< 10$.
Without any correction factor, that means that multiplier would be $1 <$ multiplier $< 10$.
When a correction factor is applied, $1.4$ in this case, the range shifts to $1.4 <$ multiplier $< 14$.

Let's follow the example that we fall slightly short of the min time, 0.9x of the min time. Without any correction factor, we could just hit the 1.0x of min time with (roughly) multiplier $= 1/0.9 \simeq 1.11 \ll 1.4$.

[LebedevRI]

According to my calculations, in step 4 you would run somewhere between 1.4 × N and 14 × N .

Since GetMinTimeToApply() > i.seconds, we always have GetMinTimeToApply() / i.seconds > 1 . Given that i.seconds / GetMinTimeToApply() > 0.1 , we have GetMinTimeToApply() / i.seconds < 10 . Without any correction factor, that means that multiplier would be 1 < multiplier < 10 . When a correction factor is applied, 1.4 in this case, the range shifts to 1.4 < multiplier < 14 .

Right.

Let's follow the example that we fall slightly short of the min time, 0.9x of the min time. Without any correction factor, we could just hit the 1.0x of min time with (roughly) multiplier = 1 / 0.9 ≃ 1.11 ≪ 1.4 .

But you do see the problem, right? We will then have spent 1.9x min_time.

[LebedevRI]

Very quick-n-dirty back of the napkin proof:

julia code

using Unitful
using Distributions
using Random

MIN_TIME = 1u"s"
IDEAL_NUM_ITERATIONS = 1000
ITERATION_MEAN_TIME = MIN_TIME / IDEAL_NUM_ITERATIONS

ITERATION_CoV = 0.01
ITERATION_STDEV = ITERATION_CoV * ITERATION_MEAN_TIME

ITERATION_STEP = 10
# TIME_FUDJE = 1.1

SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = 0.25
SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = 0.10

function run_test(;TIME_FUDJE) 
    iterations_total = missing
    time_total = missing
    
    iterations = missing
    time = missing

    while true
        prev_iterations = iterations
        if iterations isa Missing
            iterations = 1
        elseif (time / MIN_TIME) <= (ITERATION_STEP/100)
            iterations *= ITERATION_STEP
        else
            multiplier = MIN_TIME * TIME_FUDJE / max(time, 1e-9u"s")
            @show multiplier
            iterations *= multiplier
        end
        iterations = convert(Int64, ceil(iterations))
        @assert (prev_iterations isa Missing) || iterations > prev_iterations
        
        ds = Bernoulli(SYSTEM_BACKGROUND_JITTER_LIKELYHOOD)
        di = Normal(ustrip(u"s", ITERATION_MEAN_TIME), ustrip(u"s", ITERATION_STDEV))
        
        background_jitter_scale = [ (!b ? 1 : (1 + SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)) for b in rand(ds, iterations)]
        iteration_times = background_jitter_scale .* (rand(di, iterations).*u"s")
        time = sum(iteration_times)
        
        time_total = sum(skipmissing([time, time_total]))
        iterations_total = sum(skipmissing([iterations, iterations_total]))

        if time > MIN_TIME
            break
        end
    end
    @assert time > MIN_TIME
    return (iterations_total, time_total)
end

@show run_test(TIME_FUDJE=1.4)
@show run_test(TIME_FUDJE=1.0)

multiplier = 13.672021839106499
run_test(TIME_FUDJE = 1.4) = (1479, 1.5175746366315466 s)
multiplier = 9.686340899121179
multiplier = 1.0067496719348774
run_test(TIME_FUDJE = 1.0) = (2056, 2.1082144456986898 s)

So 1.0 fudge factor ends up being slower.

NOTE: that model is not meant to be authoritative.
I'm not sure if modelling system background load as a Bernoulli distribution
is the right approach, but i do believe that it is a right first-order approximation.

[LebedevRI]

Some plotting:

julia code

module MyUnits
    using Unitful
    #export pt, mm, cm
    @unit fudgeFactor "×" FudgeFactor 1 false
    @unit iterations "iters" Iterations 1 false

    #const mm = u"mm"
    #const cm = u"cm"

    Unitful.register(@__MODULE__)
    function __init__()
        return Unitful.register(@__MODULE__)
    end               
end

using Unitful
using Distributions
using Random
using Measurements
using Plots

plotlyjs()

MIN_TIME = 1u"s"
IDEAL_NUM_ITERATIONS = 100
ITERATION_MEAN_TIME = MIN_TIME / IDEAL_NUM_ITERATIONS

ITERATION_CoV = 1u"percent"
ITERATION_STDEV = ITERATION_CoV * ITERATION_MEAN_TIME

ITERATION_STEP = 10
AUTHORITATIVE_TIME = 10u"percent"

SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = 25u"percent"
SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = 0.10

NUM_REPETITIONS = 1000

function run_test_repetition(;TIME_FUDJE) 
    iterations_total = missing
    time_total = missing
    
    iterations = missing
    time = missing

    while true
        prev_iterations = iterations
        if iterations isa Missing
            iterations = 1
        elseif ((time / MIN_TIME)*100u"percent") <= AUTHORITATIVE_TIME
            iterations *= ITERATION_STEP
        else
            multiplier = MIN_TIME * TIME_FUDJE / max(time, 1e-9u"s")
            iterations *= multiplier
        end
        iterations = convert(Int64, ceil(iterations))
        @assert (prev_iterations isa Missing) || iterations > prev_iterations
        
        ds = Bernoulli(ustrip(u"percent", SYSTEM_BACKGROUND_JITTER_LIKELYHOOD)/100)
        di = Normal(ustrip(u"s", ITERATION_MEAN_TIME), ustrip(u"s", ITERATION_STDEV))
        
        background_jitter_scale = [ (!b ? 1 : (1 + SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)) for b in rand(ds, iterations)]
        iteration_times = background_jitter_scale .* (rand(di, iterations).*u"s")
        time = sum(iteration_times)
        
        iterations_total = sum(skipmissing([iterations, iterations_total]))
        time_total = sum(skipmissing([time, time_total]))

        if time > MIN_TIME
            break
        end
    end
    @assert time > MIN_TIME
    return (iterations_total, time_total)
end

function run_test(TIME_FUDJE)
    v = [run_test_repetition(TIME_FUDJE=TIME_FUDJE) for r in 1:NUM_REPETITIONS]
    all_iterations = getfield.(v, 1)
    all_times = getfield.(v, 2)
   
    iterations_mean = mean(all_iterations)
    time_mean = mean(all_times)
    
    iterations_std = std(all_iterations; mean=iterations_mean)
    time_std = std(all_times; mean=time_mean)
    
    return (measurement(iterations_mean, iterations_std),
            measurement(time_mean, time_std))
end

FUDGE_MIN = 1.0
FUDGE_MAX = 1.5
FUDGE_STEP = 0.01

TIME_FUDJES = collect(range(FUDGE_MIN,stop=FUDGE_MAX,length=1+ceil(Int64, (FUDGE_MAX-FUDGE_MIN)/FUDGE_STEP)))
RESULTS = run_test.(TIME_FUDJES)
TIME_FUDJES = TIME_FUDJES.*u"fudgeFactor"
# ITERATIONS = getfield.(RESULTS, 1).*u"iterations"
TIMES = getfield.(RESULTS, 2)

plot(TIME_FUDJES, TIMES, xlabel="fudge factor", ylabel="total time", size=(800,600))

Without touching the rest of the params there ^,
with 10ms iterations (x 100 repetitions), we get:

but with 1ms iterations (x 1000 repetitions), we get:

Maybe the factor of 1.4 is too much, but i'm certain that 1.0 is just wrong.

[LebedevRI]

And some more fancy-ness (most things are now sampled from distributions
instead of being hardcoded. again, not 100% confident in particular choices):

julia code

module MyUnits
    using Unitful

    @unit fudgeFactor "×" FudgeFactor 1 false
    @unit iterations "iters" Iterations 1 false

    Unitful.register(@__MODULE__)
    function __init__()
        return Unitful.register(@__MODULE__)
    end               
end

using Base.Threads
using Base.Iterators
using Unitful
using Distributions
using Random
using Measurements
using ProgressMeter
using Plots

plotlyjs()

MIN_TIME = 1u"s"
ITERATION_STEP = 10
AUTHORITATIVE_TIME = 10u"percent"

#dist_ITERATION_MEAN_TIME = LogUniform(ustrip(u"s", 1u"μs"), ustrip(u"s", 10u"s"))
#dist_ITERATION_CoV = LogUniform(0.1/100, 10/100)
#dist_SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = LogUniform(0.1/100, 10/100)
#dist_SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = LogUniform(0.1/100, 10/100)

dist_ITERATION_MEAN_TIME = censored(LogNormal(0, 2.5), ustrip(u"s", 1u"μs"), ustrip(u"s", 1u"s"))
dist_ITERATION_CoV = censored(LogNormal(1, 1), 0.1/100, 10/100)
dist_SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = censored(LogNormal(2.5, 1), 0.1/100, 10/100)
dist_SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = censored(LogNormal(2.5, 1), 0.1/100, 10/100)

function run_test_repetition_impl(;TIME_FUDJE,
                                    ITERATION_MEAN_TIME,
                                    ITERATION_STDEV,
                                    SYSTEM_BACKGROUND_JITTER_LIKELYHOOD,
                                    SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)
    dist_SYSTEM_BACKGROUND_JITTER = Bernoulli(SYSTEM_BACKGROUND_JITTER_LIKELYHOOD)
    dist_ITERATION_TIME = Normal(ITERATION_MEAN_TIME, ITERATION_STDEV)
     
    iterations_total = missing
    time_total = missing
    
    iterations = missing
    time = missing

    while true
        prev_iterations = iterations
        if iterations isa Missing
            iterations = 1
        elseif ((time / MIN_TIME)*100u"percent") <= AUTHORITATIVE_TIME
            iterations *= ITERATION_STEP
        else
            multiplier = MIN_TIME * TIME_FUDJE / max(time, 1e-9u"s")
            iterations *= multiplier
        end
        iterations = convert(Int64, ceil(iterations))
        @assert (prev_iterations isa Missing) || iterations > prev_iterations

        background_jitter_scale = [ (!b ? 1 : (1 + SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)) for b in rand(dist_SYSTEM_BACKGROUND_JITTER, iterations)]
        iteration_times = background_jitter_scale .* (rand(dist_ITERATION_TIME, iterations).*u"s")
        time = sum(iteration_times)
        
        iterations_total = sum(skipmissing([iterations, iterations_total]))
        time_total = sum(skipmissing([time, time_total]))

        if time > MIN_TIME
            break
        end
    end
    @assert time > MIN_TIME
    real_iteration_mean_time = time_total / iterations_total
    ideal_iteration_count = ceil(Int64, MIN_TIME / real_iteration_mean_time)
    ideal_run_time = real_iteration_mean_time * ideal_iteration_count
    time_overspent = time_total / ideal_run_time
    return (iterations_total, time_total, time_overspent)
end

NUM_INNER_REPETITIONS = 100

function run_test_repetition(;TIME_FUDJE) 
    ITERATION_MEAN_TIME = rand(dist_ITERATION_MEAN_TIME)
    ITERATION_CoV = rand(dist_ITERATION_CoV)
    SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = rand(dist_SYSTEM_BACKGROUND_JITTER_LIKELYHOOD)
    SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = rand(dist_SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)

    ITERATION_STDEV = ITERATION_CoV * ITERATION_MEAN_TIME

    v = [run_test_repetition_impl(TIME_FUDJE=TIME_FUDJE,
                                  ITERATION_MEAN_TIME=ITERATION_MEAN_TIME,
                                  ITERATION_STDEV=ITERATION_STDEV,
                                  SYSTEM_BACKGROUND_JITTER_LIKELYHOOD=SYSTEM_BACKGROUND_JITTER_LIKELYHOOD,
                                  SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT=SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT)
                                for r in 1:NUM_INNER_REPETITIONS]
    next!(p)
    return v
end

NUM_REPETITIONS = 10^5

function run_test(TIME_FUDJE)
    chunk_size = max(1, ceil(Int64, NUM_REPETITIONS // nthreads()))
    data_chunks = partition(1:NUM_REPETITIONS, chunk_size)
    tasks = map(data_chunks) do chunk
        @spawn begin
            state = []
            for x in chunk
                push!(state, run_test_repetition(TIME_FUDJE=TIME_FUDJE))
            end
            return state
        end
    end
    
    v = fetch.(tasks)
    v = vcat(v...)
    v = vcat(v...)
    all_iterations = getfield.(v, 1)
    all_times = getfield.(v, 2)
    all_times_overspent = getfield.(v, 3)
   
    iterations_mean = mean(all_iterations)
    time_mean = mean(all_times)
    time_overspent_mean = mean(all_times_overspent)

    iterations_std = std(all_iterations; mean=iterations_mean)
    time_std = std(all_times; mean=time_mean)
    time_overspent_std = std(all_times_overspent; mean=time_overspent_mean)
    
    return (measurement(iterations_mean, iterations_std),
            measurement(time_mean, time_std),
            measurement(time_overspent_mean, time_overspent_std))
end

FUDGE_MIN = 1.0
FUDGE_MAX = 1.5
FUDGE_STEP = 0.01

TIME_FUDJES = collect(range(FUDGE_MIN,stop=FUDGE_MAX,length=1+ceil(Int64, (FUDGE_MAX-FUDGE_MIN)/FUDGE_STEP)))

chunk_size = max(1, ceil(Int64, length(TIME_FUDJES) // nthreads()))
data_chunks = partition(TIME_FUDJES, chunk_size)
p = Progress(NUM_REPETITIONS*length(TIME_FUDJES); dt=1)
tasks = map(data_chunks) do chunk
    @spawn begin
        state = []
        for x in chunk
            push!(state, run_test(x))
        end
        return state
    end
end
RESULTS = fetch.(tasks)
RESULTS = vcat(RESULTS...)
TIME_FUDJES = TIME_FUDJES.*u"fudgeFactor"
# ITERATIONS = getfield.(RESULTS, 1).*u"iterations"
TIMES = getfield.(RESULTS, 2)
TIMES_OVERSPENT = getfield.(RESULTS, 3)

scatter(TIME_FUDJES, Measurements.value.(TIMES_OVERSPENT), #series_annotations = ann,
        xticks=TIME_FUDJES, yticks=0:0.005:3, xlabel="fudge factor", ylabel="time overspenditure (factor)", size=(1920,1080))

EDIT: the distributions there ^ are wrong, these might be better: (graphs not updated)

using SpecialFunctions

function lognormal_σ(;μ,v,num_σ=3)
    p = erf(num_σ/sqrt(2))
    return -(μ - log(v))/(sqrt(2)*erfinv(2*p-1))
end

# 0σ = 10ms, +3σ = 1s
# NOTE: in reality, mean should be 1ms, but that makes this simulation much slower...
μ=log(10/1000)
dist_ITERATION_MEAN_TIME = censored(LogNormal(μ, lognormal_σ(μ=μ,v=1000/1000)), ustrip(u"s", 1u"μs"), ustrip(u"s", 1u"s"))

# 0σ = 1%, +3σ = 25%
μ=log(1/100)
dist_ITERATION_CoV = censored(LogNormal(μ, lognormal_σ(μ=μ,v=25/100)), 0.1/100, 25/100)

# 0σ = 10%, +3σ = 25%
μ=log(10/100)
dist_SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = censored(LogNormal(μ, lognormal_σ(μ=μ,v=25/100)), 0.1/100, 50/100)

# 0σ = 5%, +3σ = 25%
μ=log(5/100)
dist_SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = censored(LogNormal(μ, lognormal_σ(μ=μ,v=25/100)), 0.1/100, 50/100)

Again, i'm not quite confident that my "statistical model" is sufficiently precise for anything,
but maybe 1.4 could be lowered somewhat, perhaps to 1.10, theoretically.

[ruben-laso]

I agree that the value 1.0 falls into the theoretical (non-real) world.

Looking at your plots, I also agree that something in the 1.02 - 1.10 range might be more sensible.

In any case, I think the best option could be to have an environment variable (or a program argument like --benchmark_min_time_thres) so the user can tweak this if needed instead of a hardcoded constant value.

[LebedevRI]

@dmah42 question: does google retain benchmark .json data for long-term performance tracking?
If it does, we maybe could at least derive (with-noise) dist_ITERATION_MEAN_TIME from that real-world data.

We don't encode not nearly enough info into that json (we are missing min_time and use_real_time),
so it's would to be somewhat lossy, but i think we could side-step that issue.
I think, only the following parts of said json output would be necessary:

{
"run_name": <?>,
"run_type": <?>,
"threads": <?>, % we were scaling by this factor until recently
"iterations": <?>,
"aggregate_name": <?>,
"real_time": <?>,
"cpu_time": <?>,
"time_unit": <?>
}

I'm guessing you'll want to censor "run_name", i think it would be fine to drop everything but /min_time:<>and/use_real_time` from that string.

Out of all the benchmarks that are tracked, have all of them been run at least once since #1836?
I think it would be least confusing to only look at the benchmarks //with// that change,
(or without it).

Is that possible?

[LebedevRI]

@dmah42 (consider my previous question to be a theoretical question, not a request.)

@ruben-laso

I agree that the value 1.0 falls into the theoretical (non-real) world.

Looking at your plots, I also agree that something in the 1.02 - 1.10 range might be more sensible.

In any case, I think the best option could be to have an environment variable (or a program argument like --benchmark_min_time_thres) so the user can tweak this if needed instead of a hardcoded constant value.

That feels very much placebo-y, as far as i can tell, even with that fudge at 1.1, we still spend 1.5x the time
we ideally need to (as per the final iteration time), which is not significantly better than what we'd spent at 1.4:

julia code

module MyUnits
    using Unitful

    @unit fudgeFactor "×" FudgeFactor 1 false
    @unit iterations "iters" Iterations 1 false

    Unitful.register(@__MODULE__)
    function __init__()
        return Unitful.register(@__MODULE__)
    end               
end

using Base.Threads
using Base.Iterators
using Unitful
using Distributions
using Random
using Measurements
using ProgressMeter
using Plots
using SpecialFunctions
using DataFrames

plotlyjs()

MIN_TIME = 1u"s"

#dist_ITERATION_MEAN_TIME = LogUniform(ustrip(u"s", 1u"μs"), ustrip(u"s", 10u"s"))
#dist_ITERATION_CoV = LogUniform(0.1/100, 10/100)
#dist_SYSTEM_BACKGROUND_JITTER_LIKELYHOOD = LogUniform(0.1/100, 10/100)
#dist_SYSTEM_BACKGROUND_JITTER_RUNTIME_IMPACT = LogUniform(0.1/100, 10/100)

function lognormal_σ(;μ,v,num_σ=3)
    p = erf(num_σ/sqrt(2))
    return -(μ - log(v))/(sqrt(2)*erfinv(2*p-1))
end

# 0σ = 1ms, +3σ = 1s
# NOTE: in reality, mean should be 1ms, but that makes this simulation much slower...
μ=log(1/1000)
dist_ITERATION_MEAN_TIME = censored(LogNormal(μ, lognormal_σ(μ=μ,v=1000/1000)), ustrip(u"s", 1u"μs"), ustrip(u"s", 1u"s"))

# 0σ = 10%, +3σ = 25%
μ=log(10/100)
dist_ITERATION_CoV = censored(LogNormal(μ, lognormal_σ(μ=μ,v=25/100)), 0.1/100, 25/100)

struct BenchmarkResult
    AUTHORITATIVE_TIME::Float64
    ITERATION_COUNT_MULTIPLIER::Int64
    TIME_FUDJE::Float64
    time_overspent::Any
end

function run_test_repetition_impl(; AUTHORITATIVE_TIME,
                                    ITERATION_COUNT_MULTIPLIER,
                                    TIME_FUDJE)
    ITERATION_MEAN_TIME = rand(dist_ITERATION_MEAN_TIME)
    ITERATION_CoV = rand(dist_ITERATION_CoV)

    ITERATION_STDEV = ITERATION_CoV * ITERATION_MEAN_TIME

    dist_ITERATION_TIME = Normal(ITERATION_MEAN_TIME, ITERATION_STDEV)
     
    iterations_total = missing
    time_total = missing
    
    iterations = missing
    time = missing

    while true
        prev_iterations = iterations
        if iterations isa Missing
            iterations = 1
        elseif (time / MIN_TIME) <= AUTHORITATIVE_TIME
            iterations *= ITERATION_COUNT_MULTIPLIER
        else
            multiplier = MIN_TIME * TIME_FUDJE / max(time, 1e-9u"s")
            iterations *= multiplier
        end
        iterations = convert(Int64, ceil(iterations))
        @assert (prev_iterations isa Missing) || iterations > prev_iterations

        iteration_times = rand(dist_ITERATION_TIME,iterations).*u"s"
        time = sum(iteration_times)
        
        iterations_total = sum(skipmissing([iterations, iterations_total]))
        time_total = sum(skipmissing([time, time_total]))

        if time >= MIN_TIME
            break
        end
    end
    @assert time > MIN_TIME
    real_iteration_mean_time = time_total / iterations_total
    ideal_iteration_count = ceil(Int64, MIN_TIME / real_iteration_mean_time)
    ideal_run_time = real_iteration_mean_time * ideal_iteration_count
    time_overspent = time_total / ideal_run_time
    return time_overspent
end

NUM_BENCHMARK_REPETITIONS = 10^6

function run_test_configuration(; AUTHORITATIVE_TIME,
                                  ITERATION_COUNT_MULTIPLIER,
                                  TIME_FUDJE) 
    chunk_size = max(1, ceil(Int64, NUM_BENCHMARK_REPETITIONS // nthreads()))
    data_chunks = partition(1:NUM_BENCHMARK_REPETITIONS, chunk_size)
    tasks = map(data_chunks) do chunk
        @spawn begin
            state = []
            for x in chunk
                push!(state, run_test_repetition_impl(
                                  AUTHORITATIVE_TIME=AUTHORITATIVE_TIME,
                                  ITERATION_COUNT_MULTIPLIER=ITERATION_COUNT_MULTIPLIER,
                                  TIME_FUDJE=TIME_FUDJE))
                next!(p)
            end
            return state
        end
    end
    
    v = fetch.(tasks)
    v = vcat(v...)
    v = vcat(v...)

    time_overspent = v

    time_overspent_mean = mean(time_overspent)
    time_overspent_std = std(time_overspent; mean=time_overspent_mean)
    time_overspent = measurement(time_overspent_mean, time_overspent_std)
    #time_overspent = time_overspent_mean + 3*time_overspent_std
    
    v = BenchmarkResult(AUTHORITATIVE_TIME,
                        ITERATION_COUNT_MULTIPLIER,
                        TIME_FUDJE,
                        time_overspent)
    return v
end

function run_test(CONFIGURATIONS)    
    chunk_size = max(1, ceil(Int64, length(CONFIGURATIONS) // nthreads()))
    data_chunks = partition(CONFIGURATIONS, chunk_size)
    tasks = map(data_chunks) do chunk
        @spawn begin
            state = []
            for x in chunk
                push!(state, run_test_configuration(AUTHORITATIVE_TIME=x[1],ITERATION_COUNT_MULTIPLIER=x[2],TIME_FUDJE=x[3]))
            end
            return state
        end
    end

    v = fetch.(tasks)
    v = vcat(v...)
    v = vcat(v...)
    return v
end

#AUTHORITATIVE_TIME = [10/100]
#ITERATION_COUNT_MULTIPLIER = [10]
#TIME_FUDGE = [1.4]

AUTHORITATIVE_TIME = [10/100]
ITERATION_COUNT_MULTIPLIER = [10]
TIME_FUDGE = [1.0,1.1,1.4]

#AUTHORITATIVE_TIME = 1/100:1/100:5/100
#ITERATION_COUNT_MULTIPLIER = 2:1:20
#TIME_FUDGE = 1.0:0.01:1.1

CONFIGURATIONS = [(x,y,z) for x in AUTHORITATIVE_TIME
                          for y in ITERATION_COUNT_MULTIPLIER
                          for z in TIME_FUDGE]
NUM_CONFIGURATIONS = length(CONFIGURATIONS)

p = Progress(NUM_BENCHMARK_REPETITIONS*NUM_CONFIGURATIONS; dt=1)
R = run_test(CONFIGURATIONS)

X = [r.ITERATION_COUNT_MULTIPLIER for r in R]
Y = [r.AUTHORITATIVE_TIME for r in R]
Y2 = [r.TIME_FUDJE for r in R]
Z = [r.time_overspent for r in R]

df = DataFrame(ITERATION_COUNT_MULTIPLIER=X, AUTHORITATIVE_TIME=Y, TIME_FUDJE=Y2, time_overspent=Z)
sort!(df, [:time_overspent])
sort!(df, [order(:time_overspent), order(:AUTHORITATIVE_TIME), order(:ITERATION_COUNT_MULTIPLIER, rev=true), order(:TIME_FUDJE, rev=true)])
@show first(df, 100)
@assert false

Z = Measurements.value.(Z)

#heatmap(X,Y,Z; size=(800,600),
#        xlabel="ITERATION_COUNT_MULTIPLIER (x)", ylabel="AUTHORITATIVE_TIME (%)", zlabel="time_overspent (x)",
#        xticks=ITERATION_COUNT_MULTIPLIER, yticks=AUTHORITATIVE_TIME)

scatter(X,Y,Z; zcolor=Z, size=(800,600),
        xlabel="ITERATION_COUNT_MULTIPLIER (x)", ylabel="AUTHORITATIVE_TIME (%)", zlabel="time_overspent (x)",
        xticks=ITERATION_COUNT_MULTIPLIER, yticks=AUTHORITATIVE_TIME)

#scatter(Y2,Z; size=(800,600),
#        xlabel="TIME_FUDGE", ylabel="time_overspent (x)",
#        xticks=TIME_FUDGE)

first(df, 100) = 3×4 DataFrame
 Row │ ITERATION_COUNT_MULTIPLIER  AUTHORITATIVE_TIME  TIME_FUDJE  time_overspent
     │ Int64                       Float64             Float64     Measurement…
─────┼────────────────────────────────────────────────────────────────────────────
   1 │                         10                 0.1         1.1       1.53±0.28
   2 │                         10                 0.1         1.4       1.82±0.28
   3 │                         10                 0.1         1.0       2.04±0.85

We could also lower AUTHORITATIVE_TIME and maybe change ITERATION_COUNT_MULTIPLIER,
and lower the total time spent to as low as 1.2x of the ideal,
but really does feel ad-hoc to change anything without better statistical model...

[dmah42]

the right solution to this is to run until we reduce the confidence interval for the iterations within a run to a particular threshold (or max time to avoid infinite runs). anything other than this is just fine-tuning and is unlikely to scale across all the various use-cases.

if someone wants to reduce 1.4 to 1.2 or 1.1 or whatever then create a PR but i really don't think it matters all that much (and thank you @LebedevRI for basically proving that :D ).

[LebedevRI]

the right solution to this is to run until we reduce the confidence interval for the iterations within a run to a particular threshold (or max time to avoid infinite runs). anything other than this is just fine-tuning and is unlikely to scale across all the various use-cases.

if someone wants to reduce 1.4 to 1.2 or 1.1 or whatever then create a PR but i really don't think it matters all that much (and thank you @LebedevRI for basically proving that :D ).

Yeah, that's my general understanding, too. I haven't yet tried modelling said different convergence approach.
Trivia bit: this time_overspent factor for current implementation seems to best described by inverse gaussian distribution.