22882410591

Committed 10 Mar 2026 01:12AM UTC coverage: 94.597%. Remained the same

Build # 22882410591

Build Type

push

github

Committed by

PieterjanRobbe

Commit Message

Removed manifest file

Coverage Stats

893 of 944 relevant lines covered (94.6%)

289082.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.86

/src/methods.jl

## methods.jl : refactored methods from run.jl
#
# Implements refactored methods from run.jl.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods

#
# basics
#
ndims(::Estimator{<:AbstractIndexSet{d}}) where d = d

get_index_set(estimator::Estimator, sz) = get_index_set(estimator.index_set, sz)

get_tols(estimator::Estimator, tol::T) where T<:Real = estimator[:continuate] ? estimator[:continuation_mul_factor].^(estimator[:nb_of_tols]-1:-1:0)*tol : T[tol]

mse(estimator::Estimator) = varest(estimator) + bias(estimator)^2

rmse(estimator::Estimator) = sqrt(mse(estimator))

converged(estimator::Estimator, ϵ::Real, θ::Real) = ( bias(estimator)^2 ≤ (1-θ)*ϵ^2 || mse(estimator) ≤ ϵ^2 )

max_level_exceeded(estimator::Estimator) = sz(estimator) ≥ estimator[:max_index_set_param]

#
# inspector functions: mean, var, varest...
#
function aggregate(estimator::Estimator, values)
    if estimator[:aggregation] == "sum"
        return sum(values)
    elseif estimator[:aggregation] == "norm"
        return norm(values)
    else
        if estimator[:qoi_with_max_var] == 0
            qoi_vars = map(qoi -> var(estimator, qoi), 1:estimator[:nb_of_qoi])
            qoi_with_max_var = argmax(qoi_vars)
        else
            qoi_with_max_var = estimator[:qoi_with_max_var]
        end
        return values[qoi_with_max_var]
    end
end

var(estimator::Estimator{<:AbstractIndexSet, <:MC}, qoi::Integer) = sum(var(samples_diff(estimator, qoi, index)) for index in keys(estimator))

var(estimator::Estimator{<:AbstractIndexSet, <:QMC}, qoi::Integer) = sum(mean(var(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)) for index in keys(estimator))

for f in [:mean, :var]
    @eval begin

    $f(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = aggregate(estimator, map(qoi -> $f(samples_diff(estimator, qoi, index)), 1:estimator[:nb_of_qoi]))

    $(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = aggregate(estimator, map(qoi -> $f(samples(estimator, qoi, index)), 1:estimator[:nb_of_qoi]))

    $f(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> mean($f(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)), 1:estimator[:nb_of_qoi]))

    $(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> mean($f(samples(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)), 1:estimator[:nb_of_qoi]))

    end
end

mean(estimator::Estimator) = sum(mean(estimator, index) for index in keys(estimator))

cost(estimator::Estimator, index::Index) = estimator[:cost_model] isa EmptyFunction ? time(estimator, index) : work(estimator, index)

varest(estimator::Estimator, index::Index) = var(estimator, index) / nb_of_samples(estimator, index)

varest(estimator::Estimator) =  sum(varest(estimator, index) for index in keys(estimator))

#
# rates
#
for (f, g, sgn) in zip([:α, :β, :γ], [:mean, :var, :cost], [-1, -1, 1])
    @eval begin

        $f(estimator::Estimator{<:SL}) = nothing

        $f(estimator::Estimator{<:AbstractIndexSet}) = $sgn.*getindex.(broadcast(i->$(Symbol("rates_", f))(estimator, i), 1:ndims(estimator)), 2)

        $(Symbol("rates_", f))(estimator::Estimator{<:AbstractML}) = $(Symbol("rates_", f))(estimator, 1)

        $(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, dir::Integer) = $(Symbol("rates_", f))(estimator, (maximum(getindex.(keys(estimator), dir)) + 1) * Index(ntuple(i -> i == dir, ndims(estimator))), dir)

        function $(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, idx::Index, dir::Integer)
            m = idx[dir] - 1
            if m < 2
                return (NaN, NaN)
            else
                x = m:-1:1
                y = map(xᵢ -> $g(estimator, idx - xᵢ * Index(ntuple(i -> i == dir, ndims(estimator)))), 1:m)
                idcs = .!(isnan.(y) .| iszero.(y))
                θ = interp1(view(x, idcs), log2.(abs.(view(y, idcs))))
                return tuple(θ...)
            end
        end

        function $(Symbol("_regress_", g))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index)
            p = broadcast(dir->$(Symbol("rates_", f))(estimator, index, dir), 1:ndims(estimator))
            estimates = broadcast(dir->2^(p[dir][1]+index[dir]*p[dir][2]), 1:ndims(estimator))
            estimate = mean(filter(!isnan, estimates))
            if isnan(estimate)
                p = interp($g, estimator)
                return 2^(p[1]+sum(p[2:end].*index.I))
            else
                return estimate
            end
        end
    end
end

#
# regression
#
function interp1(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    A = hcat(ones(eltype(x), length(x)), x)
    A\y
end

function interp(f::Function, estimator::Estimator)
        filter_fcn = estimator isa Estimator{<:AbstractIndexSet, QMC} ? i -> !isempty(samples(estimator)[1][i]) : i -> length(samples(estimator)[1][i]) > 1
    idx_set = filter(filter_fcn, CartesianIndices(size(samples(estimator)[1])))
    A = [i == 0 ? 1 : getindex(index - one(index), i) for index in idx_set, i in 0:ndims(estimator)]
    y = map(i -> log2(f(estimator, i - one(i))), idx_set)
    valid = .!isinf.(y)
    A = A[valid, :]
    y = y[valid]
    try
        return A \ y
    catch e
        return fill(NaN, length(y))
    end
end

regress_mean(estimator, index) = _regress_mean(estimator, index)
regress_var(estimator, index) = _regress_var(estimator, index)
regress_cost(estimator, index) = estimator[:cost_model] isa EmptyFunction ? _regress_cost(estimator, index) : estimator[:cost_model](index)

regress_nb_of_samples(estimator::Estimator{<:SL}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])

regress_nb_of_samples(estimator::Estimator{<:SL, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])

regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)

function regress_nb_of_samples(estimator::Estimator, index_set, ϵ::Real, θ::Real, L::Integer)
    if estimator[:do_regression] && L > 2
        return _regress_nb_of_samples(estimator, index_set, ϵ, θ)
    else
        return Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)
    end
end

function _regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:MC}, index_set, ϵ::Real, θ::Real)
    vars = Dict(index => regress_var(estimator, index) for index in index_set)
    costs = Dict(index => regress_cost(estimator, index) for index in index_set)
    Σ_estimate = Σ(estimator)
    for index in keys(vars)
        Σ_estimate += sqrt(vars[index] * costs[index])
    end
    Dict(index => begin
             n_opt = optimal_nb_of_samples(ϵ, θ, vars[index], costs[index], Σ_estimate)
             max(2, min(n_opt, estimator[:nb_of_warm_up_samples]))
         end for index in index_set)
end

function compute_splitting(estimator::Estimator, ϵ::Real)
    bias_estimate = bias(estimator, max_sz(estimator))
    a = estimator[:min_splitting]
    b = estimator[:max_splitting]

    splitting = 1 - bias_estimate^2/ϵ^2
    isnan(splitting) ? a : min(b, max(a, splitting))
end

Σ(estimator::Estimator) = sum(sqrt.(map(index -> var(estimator, index) * cost(estimator, index), keys(estimator))))

optimal_nb_of_samples(estimator::Estimator, index::Index, ϵ::Real, θ::Real) = optimal_nb_of_samples(ϵ, θ, var(estimator, index), cost(estimator, index), Σ(estimator))

optimal_nb_of_samples(ϵ::Real, θ::Real, var_estimate::Real, cost_estimate::Real, Σ_estimate::Real) = ceil(Int, 1/(θ*ϵ^2) * sqrt(var_estimate/cost_estimate) * Σ_estimate)

#
# bias computation
#
boundary(estimator::Estimator, cntr::Integer) = setdiff(get_index_set(estimator, cntr), get_index_set(estimator, cntr-1))

new_index_set(estimator::Estimator, cntr::Integer) = boundary(estimator, cntr)

bias(estimator::Estimator{<:SL}) = 0.0

bias(estimator::Estimator) = bias(estimator, sz(estimator))

function bias(estimator::Estimator, sz::Integer)
    if !isempty(boundary(estimator, sz + 1) ∩ keys(estimator)) && !estimator[:robustify_bias_estimate]
        return abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, sz + 1))))
    else
        x = max(1, sz - 2):sz
        y = Float64[log2(abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, xᵢ))))) for xᵢ in x]
        idcs = isfinite.(y)
        p = interp1(x[idcs], y[idcs])
        return 2^(p[1]+(sz+1)*p[2])
    end
end

#
# adaptivity
#
profit(estimator::Estimator{<:AD}, index::Index) = abs(mean(estimator, index)) / sqrt(var(estimator, index) * cost(estimator, index))^estimator[:penalization]

max_level_exceeded(estimator::Estimator{<:AD}) = isempty(setdiff(get_index_set(estimator[:max_search_space], estimator[:max_index_set_param]), keys(estimator)))

function find_index_with_max_profit(estimator::Estimator{<:AD})
    indices = collect(active_set(estimator))
    profits = [profit(estimator, index) for index in indices]
    (max_profit, idx) = findmax(profits)
    if rand() < estimator[:acceptance_rate] || length(profits)==1
        max_index = indices[idx]
    else
        max_index = indices[rand(deleteat!(collect(1:length(profits)), idx))]
    end
    estimator[:verbose] && print_largest_profit(estimator, max_index, max_profit, indices, profits)
    return max_index
end

function new_index_set(estimator::Estimator{<:AD}, sz::Integer)
    d = ndims(estimator)
    if isempty(active_set(estimator))
        new_indices = Set{Index{d}}()
        max_index = Index(ntuple(i -> 0, d))
        add_to_active_set(estimator, max_index)
        push!(new_indices, max_index)
    else
        max_index = find_index_with_max_profit(estimator)
        add_to_old_set(estimator, max_index)
        remove_from_active_set(estimator, max_index)
        new_indices = Set{Index{d}}()
        push!(new_indices, max_index)
        for k in 1:d
            new_index = max_index + Index(ntuple(i -> i == k, d))
            if is_admissable(estimator, new_index)
                if new_index ∈ get_index_set(estimator[:max_search_space], estimator[:max_index_set_param])
                    add_to_active_set(estimator, new_index)
                    push!(new_indices, new_index)
                else
                    warn_max_index(estimator, new_index)
                    add_to_max_index_set(estimator, new_index)
                end
            end
        end
    end
    log_adaptive_index_set(estimator, max_index)
    return new_indices
end

function bias(estimator::Estimator{<:AD}, sz::Integer)
    b(indices) = abs(sum(broadcast(i -> mean(estimator, i), collect(indices))))
    if sz != max_sz(estimator)
        return b(active_set(estimator))
    else
        return min(b(active_set(estimator)), b(boundary(estimator)))
    end
end

#
# QMC related functions
#
function next_number_of_samples(estimator, index)
    if estimator[:sample_mul_factor] == 2
        Dict(index => nextpow(2, nb_of_samples(estimator, index) + 1))
    elseif estimator[:sample_mul_factor] ≤ 1
        Dict(index => nb_of_samples(estimator, index) + 1)
    else
        Dict(index => ceil(Int, nb_of_samples(estimator, index) * estimator[:sample_mul_factor]))
    end
end

function find_index_with_max_var_over_cost(estimator::Estimator{<:AbstractIndexSet, <:QMC})
    indices = collect(keys(estimator))
    vars = [varest(estimator, index) / cost(estimator, index) for index in indices]
    (max_var, idx) = findmax(vars)
    indices[idx]
end

varest(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> var(mean(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index), corrected=true) / estimator[:nb_of_shifts](index), 1:estimator[:nb_of_qoi]))

#
# Unbiased estimation
#
function next_number_of_samples(estimator::Estimator{<:U})
    n_total = sum(values(nb_of_samples(estimator)))
    n0 = max(estimator[:nb_of_warm_up_samples], n_total)
    n = estimator[:sample_mul_factor] ≤ 1 ? 1 : ceil(Int, n0 * (estimator[:sample_mul_factor] - 1))
    sample(pmf(estimator), n)
end

function sample(pmf::Dict{Index{d}, Float64}, n::Integer) where d
    u = rand(n)
    x = Vector{keytype(pmf)}(undef, n)
    psum = 0.
    for (index, p) in pmf
        for i in 1:length(x)
            if psum ≤ u[i] ≤ psum + p
                x[i] = index
            end
        end
        psum = psum + p
    end
    Dict(index => sum(x .== Ref(index)) for index in keys(pmf))
end

Geometric(p, k) = (1 - p)^k*p

function update_pmf(estimator::Estimator{<:U})
    f(estimator, index) = sqrt(var(estimator, index)/cost(estimator, index))
    p = interp(f, estimator)
    if !any(isnan.(p)) && all(p[2:end] .< 0)
        for index in keys(estimator)
            if isnan(f(estimator, index)) || isinf(f(estimator, index))
                set_pmf_key(estimator, index, 2^(p[1]+sum(p[2:end].*index.I)))
            else
                set_pmf_key(estimator, index, f(estimator, index))
            end
        end
        normalize!(pmf(estimator))
    end
end

function normalize!(pmf::Dict{<:Index, Float64})
    tot_sum = sum(values(pmf))
    for (key, val) in pmf
        pmf[key] /= tot_sum
    end
end

Prob(estimator::Estimator{<:U}, index::Index) = sum(val for (key, val) in pmf(estimator) if key ≥ index)

Prob(estimator::Estimator{<:U}) = Dict(index => Prob(estimator, index) for index in keys(estimator))

var(estimator::Estimator{<:U, <:MC}, qoi::Integer) = var(accumulator(estimator, qoi))

var(estimator::Estimator{<:U, <:QMC}, qoi::Integer) = mean(var(accumulator(estimator, qoi, shift)) for shift in 1:n_shifts(estimator))

varest(estimator::Estimator{<:U, <:MC}) = aggregate(estimator, map(qoi -> var(accumulator(estimator, qoi)) / length(accumulator(estimator, qoi)), 1:estimator[:nb_of_qoi]))

varest(estimator::Estimator{<:U, <:QMC}) = aggregate(estimator, map(qoi -> var(mean(accumulator(estimator, qoi, shift)) for shift in 1:n_shifts(estimator)) / n_shifts(estimator), 1:estimator[:nb_of_qoi]))

n_shifts(estimator::Estimator{<:U, <:QMC}) = size(accumulator(estimator), 2)

bias(estimator::Estimator{<:U}) = 0.0

converged(estimator::Estimator{<:U}, ϵ::Real) = mse(estimator) ≤ ϵ^2

mean(estimator::Estimator{<:U}) = aggregate(estimator, map(qoi -> mean(accumulator(estimator, qoi)), 1:estimator[:nb_of_qoi]))

1	## methods.jl : refactored methods from run.jl
2	#
3	# Implements refactored methods from run.jl.
4	#
5	# This file is part of MultilevelEstimators.jl - A Julia toolbox for
6	# Multilevel Monte Carlo Methods
7
8	#
9	# basics
10	#
11	ndims(::Estimator{<:AbstractIndexSet{d}}) where d = d	307,417✔
12
13	get_index_set(estimator::Estimator, sz) = get_index_set(estimator.index_set, sz)	325,264✔
14
15	get_tols(estimator::Estimator, tol::T) where T<:Real = estimator[:continuate] ? estimator[:continuation_mul_factor].^(estimator[:nb_of_tols]-1:-1:0)*tol : T[tol]	432✔
16
17	mse(estimator::Estimator) = varest(estimator) + bias(estimator)^2	14,700✔
18
19	rmse(estimator::Estimator) = sqrt(mse(estimator))	8,761✔
20
21	converged(estimator::Estimator, ϵ::Real, θ::Real) = ( bias(estimator)^2 ≤ (1-θ)*ϵ^2 \|\| mse(estimator) ≤ ϵ^2 )	12,968✔
22
23	max_level_exceeded(estimator::Estimator) = sz(estimator) ≥ estimator[:max_index_set_param]	8,168✔
24
25	#
26	# inspector functions: mean, var, varest...
27	#
28	function aggregate(estimator::Estimator, values)	1,797,335✔
29	if estimator[:aggregation] == "sum"	1,797,335✔
30	return sum(values)	595,757✔
31	elseif estimator[:aggregation] == "norm"	1,201,578✔
32	return norm(values)	585,757✔
33	else
34	if estimator[:qoi_with_max_var] == 0	615,821✔
35	qoi_vars = map(qoi -> var(estimator, qoi), 1:estimator[:nb_of_qoi])	×
36	qoi_with_max_var = argmax(qoi_vars)	×
37	else
38	qoi_with_max_var = estimator[:qoi_with_max_var]	615,821✔
39	end
40	return values[qoi_with_max_var]	615,821✔
41	end
42	end
43
44	var(estimator::Estimator{<:AbstractIndexSet, <:MC}, qoi::Integer) = sum(var(samples_diff(estimator, qoi, index)) for index in keys(estimator))	×
45
46	var(estimator::Estimator{<:AbstractIndexSet, <:QMC}, qoi::Integer) = sum(mean(var(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)) for index in keys(estimator))	×
47
48	for f in [:mean, :var]
49	@eval begin
50
51	$f(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = aggregate(estimator, map(qoi -> $f(samples_diff(estimator, qoi, index)), 1:estimator[:nb_of_qoi]))	18,587,774✔
52
53	$(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = aggregate(estimator, map(qoi -> $f(samples(estimator, qoi, index)), 1:estimator[:nb_of_qoi]))	3,932,884✔
54
55	$f(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> mean($f(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)), 1:estimator[:nb_of_qoi]))	6,100,962✔
56
57	$(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> mean($f(samples(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index)), 1:estimator[:nb_of_qoi]))	2,614,752✔
58
59	end
60	end
61
62	mean(estimator::Estimator) = sum(mean(estimator, index) for index in keys(estimator))	3,240✔
63
64	cost(estimator::Estimator, index::Index) = estimator[:cost_model] isa EmptyFunction ? time(estimator, index) : work(estimator, index)	522,179✔
65
66	varest(estimator::Estimator, index::Index) = var(estimator, index) / nb_of_samples(estimator, index)	92,560✔
67
68	varest(estimator::Estimator) = sum(varest(estimator, index) for index in keys(estimator))	29,730✔
69
70	#
71	# rates
72	#
73	for (f, g, sgn) in zip([:α, :β, :γ], [:mean, :var, :cost], [-1, -1, 1])
74	@eval begin
75
76	$f(estimator::Estimator{<:SL}) = nothing	×
77
78	$f(estimator::Estimator{<:AbstractIndexSet}) = $sgn.*getindex.(broadcast(i->$(Symbol("rates_", f))(estimator, i), 1:ndims(estimator)), 2)	93,483✔
79
80	$(Symbol("rates_", f))(estimator::Estimator{<:AbstractML}) = $(Symbol("rates_", f))(estimator, 1)	×
81
82	$(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, dir::Integer) = $(Symbol("rates_", f))(estimator, (maximum(getindex.(keys(estimator), dir)) + 1) * Index(ntuple(i -> i == dir, ndims(estimator))), dir)	155,430✔
83
84	function $(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, idx::Index, dir::Integer)	82,686✔
85	m = idx[dir] - 1	82,686✔
86	if m < 2	82,686✔
87	return (NaN, NaN)	17,286✔
88	else
89	x = m:-1:1	130,800✔
90	y = map(xᵢ -> $g(estimator, idx - xᵢ * Index(ntuple(i -> i == dir, ndims(estimator)))), 1:m)	502,558✔
91	idcs = .!(isnan.(y) .\| iszero.(y))	65,400✔
92	θ = interp1(view(x, idcs), log2.(abs.(view(y, idcs))))	118,054✔
93	return tuple(θ...)	65,400✔
94	end
95	end
96
97	function $(Symbol("_regress_", g))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index)	11,610✔
98	p = broadcast(dir->$(Symbol("rates_", f))(estimator, index, dir), 1:ndims(estimator))	34,290✔
99	estimates = broadcast(dir->2^(p[dir][1]+index[dir]*p[dir][2]), 1:ndims(estimator))	34,290✔
100	estimate = mean(filter(!isnan, estimates))	19,372✔
101	if isnan(estimate)	11,610✔
102	p = interp($g, estimator)	2,688✔
103	return 2^(p[1]+sum(p[2:end].*index.I))	2,688✔
104	else
105	return estimate	8,922✔
106	end
107	end
108	end
109	end
110
111	#
112	# regression
113	#
114	function interp1(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})	97,509✔
115	A = hcat(ones(eltype(x), length(x)), x)	255,571✔
116	A\y	105,256✔
117	end
118
119	function interp(f::Function, estimator::Estimator)	3,043✔
120	filter_fcn = estimator isa Estimator{<:AbstractIndexSet, QMC} ? i -> !isempty(samples(estimator)[1][i]) : i -> length(samples(estimator)[1][i]) > 1	49,859✔
121	idx_set = filter(filter_fcn, CartesianIndices(size(samples(estimator)[1])))	3,043✔
122	A = [i == 0 ? 1 : getindex(index - one(index), i) for index in idx_set, i in 0:ndims(estimator)]	3,043✔
123	y = map(i -> log2(f(estimator, i - one(i))), idx_set)	24,618✔
124	valid = .!isinf.(y)	3,043✔
125	A = A[valid, :]	3,043✔
126	y = y[valid]	3,043✔
127	try	3,043✔
128	return A \ y	3,043✔
129	catch e
130	return fill(NaN, length(y))	5✔
131	end
132	end
133
134	regress_mean(estimator, index) = _regress_mean(estimator, index)	×
135	regress_var(estimator, index) = _regress_var(estimator, index)	5,805✔
136	regress_cost(estimator, index) = estimator[:cost_model] isa EmptyFunction ? _regress_cost(estimator, index) : estimator[:cost_model](index)	5,805✔
137
138	regress_nb_of_samples(estimator::Estimator{<:SL}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])	×
139
140	regress_nb_of_samples(estimator::Estimator{<:SL, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])	×
141
142	regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)	6,482✔
143
144	function regress_nb_of_samples(estimator::Estimator, index_set, ϵ::Real, θ::Real, L::Integer)	6,482✔
145	if estimator[:do_regression] && L > 2	6,482✔
146	return _regress_nb_of_samples(estimator, index_set, ϵ, θ)	1,622✔
147	else
148	return Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)	4,860✔
149	end
150	end
151
152	function _regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:MC}, index_set, ϵ::Real, θ::Real)	1,622✔
153	vars = Dict(index => regress_var(estimator, index) for index in index_set)	1,622✔
154	costs = Dict(index => regress_cost(estimator, index) for index in index_set)	1,622✔
155	Σ_estimate = Σ(estimator)	1,622✔
156	for index in keys(vars)	1,622✔
157	Σ_estimate += sqrt(vars[index] * costs[index])	5,805✔
158	end	5,805✔
159	Dict(index => begin	1,622✔
160	n_opt = optimal_nb_of_samples(ϵ, θ, vars[index], costs[index], Σ_estimate)	5,805✔
161	max(2, min(n_opt, estimator[:nb_of_warm_up_samples]))	5,805✔
162	end for index in index_set)
163	end
164
165	function compute_splitting(estimator::Estimator, ϵ::Real)	13,652✔
166	bias_estimate = bias(estimator, max_sz(estimator))	13,652✔
167	a = estimator[:min_splitting]	13,652✔
168	b = estimator[:max_splitting]	13,652✔
169
170	splitting = 1 - bias_estimate^2/ϵ^2	13,652✔
171	isnan(splitting) ? a : min(b, max(a, splitting))	13,652✔
172	end
173
174	Σ(estimator::Estimator) = sum(sqrt.(map(index -> var(estimator, index) * cost(estimator, index), keys(estimator))))	288,168✔
175
176	optimal_nb_of_samples(estimator::Estimator, index::Index, ϵ::Real, θ::Real) = optimal_nb_of_samples(ϵ, θ, var(estimator, index), cost(estimator, index), Σ(estimator))	31,195✔
177
178	optimal_nb_of_samples(ϵ::Real, θ::Real, var_estimate::Real, cost_estimate::Real, Σ_estimate::Real) = ceil(Int, 1/(θϵ^2) sqrt(var_estimate/cost_estimate) * Σ_estimate)	37,000✔
179
180	#
181	# bias computation
182	#
183	boundary(estimator::Estimator, cntr::Integer) = setdiff(get_index_set(estimator, cntr), get_index_set(estimator, cntr-1))	162,632✔
184
185	new_index_set(estimator::Estimator, cntr::Integer) = boundary(estimator, cntr)	10,800✔
186
187	bias(estimator::Estimator{<:SL}) = 0.0	×
188
189	bias(estimator::Estimator) = bias(estimator, sz(estimator))	34,152✔
190
191	function bias(estimator::Estimator, sz::Integer)	39,856✔
192	if !isempty(boundary(estimator, sz + 1) ∩ keys(estimator)) && !estimator[:robustify_bias_estimate]	39,856✔
193	return abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, sz + 1))))	×
194	else
195	x = max(1, sz - 2):sz	40,126✔
196	y = Float64[log2(abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, xᵢ))))) for xᵢ in x]	443,479✔
197	idcs = isfinite.(y)	39,856✔
198	p = interp1(x[idcs], y[idcs])	39,856✔
199	return 2^(p[1]+(sz+1)*p[2])	39,856✔
200	end
201	end
202
203	#
204	# adaptivity
205	#
206	profit(estimator::Estimator{<:AD}, index::Index) = abs(mean(estimator, index)) / sqrt(var(estimator, index) * cost(estimator, index))^estimator[:penalization]	2,709✔
207
208	max_level_exceeded(estimator::Estimator{<:AD}) = isempty(setdiff(get_index_set(estimator[:max_search_space], estimator[:max_index_set_param]), keys(estimator)))	1,624✔
209
210	function find_index_with_max_profit(estimator::Estimator{<:AD})	1,624✔
211	indices = collect(active_set(estimator))	2,164✔
212	profits = [profit(estimator, index) for index in indices]	1,624✔
213	(max_profit, idx) = findmax(profits)	1,624✔
214	if rand() < estimator[:acceptance_rate] \|\| length(profits)==1	1,624✔
215	max_index = indices[idx]	1,624✔
216	else
217	max_index = indices[rand(deleteat!(collect(1:length(profits)), idx))]	×
218	end
219	estimator[:verbose] && print_largest_profit(estimator, max_index, max_profit, indices, profits)	1,624✔
220	return max_index	1,624✔
221	end
222
223	function new_index_set(estimator::Estimator{<:AD}, sz::Integer)	2,164✔
224	d = ndims(estimator)	2,164✔
225	if isempty(active_set(estimator))	2,164✔
226	new_indices = Set{Index{d}}()	540✔
227	max_index = Index(ntuple(i -> 0, d))	900✔
228	add_to_active_set(estimator, max_index)	540✔
229	push!(new_indices, max_index)	540✔
230	else
231	max_index = find_index_with_max_profit(estimator)	1,624✔
232	add_to_old_set(estimator, max_index)	1,624✔
233	remove_from_active_set(estimator, max_index)	1,624✔
234	new_indices = Set{Index{d}}()	1,624✔
235	push!(new_indices, max_index)	1,624✔
236	for k in 1:d	1,624✔
237	new_index = max_index + Index(ntuple(i -> i == k, d))	8,660✔
238	if is_admissable(estimator, new_index)	6,496✔
239	if new_index ∈ get_index_set(estimator[:max_search_space], estimator[:max_index_set_param])	2,297✔
240	add_to_active_set(estimator, new_index)	2,296✔
241	push!(new_indices, new_index)	2,296✔
242	else
243	warn_max_index(estimator, new_index)	1✔
244	add_to_max_index_set(estimator, new_index)	1✔
245	end
246	end
247	end	3,248✔
248	end
249	log_adaptive_index_set(estimator, max_index)	2,164✔
250	return new_indices	2,164✔
251	end
252
253	function bias(estimator::Estimator{<:AD}, sz::Integer)	7,948✔
254	b(indices) = abs(sum(broadcast(i -> mean(estimator, i), collect(indices))))	54,583✔
255	if sz != max_sz(estimator)	7,948✔
256	return b(active_set(estimator))	1,087✔
257	else
258	return min(b(active_set(estimator)), b(boundary(estimator)))	6,861✔
259	end
260	end
261
262	#
263	# QMC related functions
264	#
265	function next_number_of_samples(estimator, index)	688✔
266	if estimator[:sample_mul_factor] == 2	688✔
267	Dict(index => nextpow(2, nb_of_samples(estimator, index) + 1))	×
268	elseif estimator[:sample_mul_factor] ≤ 1	688✔
269	Dict(index => nb_of_samples(estimator, index) + 1)	×
270	else
271	Dict(index => ceil(Int, nb_of_samples(estimator, index) * estimator[:sample_mul_factor]))	688✔
272	end
273	end
274
275	function find_index_with_max_var_over_cost(estimator::Estimator{<:AbstractIndexSet, <:QMC})	688✔
276	indices = collect(keys(estimator))	1,146✔
277	vars = [varest(estimator, index) / cost(estimator, index) for index in indices]	688✔
278	(max_var, idx) = findmax(vars)	688✔
279	indices[idx]	688✔
280	end
281
282	varest(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = aggregate(estimator, map(qoi -> var(mean(samples_diff(estimator, qoi, shift, index)) for shift in 1:estimator[:nb_of_shifts](index), corrected=true) / estimator[:nb_of_shifts](index), 1:estimator[:nb_of_qoi]))	1,851,178✔
283
284	#
285	# Unbiased estimation
286	#
287	function next_number_of_samples(estimator::Estimator{<:U})	355✔
288	n_total = sum(values(nb_of_samples(estimator)))	635✔
289	n0 = max(estimator[:nb_of_warm_up_samples], n_total)	355✔
290	n = estimator[:sample_mul_factor] ≤ 1 ? 1 : ceil(Int, n0 * (estimator[:sample_mul_factor] - 1))	710✔
291	sample(pmf(estimator), n)	355✔
292	end
293
294	function sample(pmf::Dict{Index{d}, Float64}, n::Integer) where d	355✔
295	u = rand(n)	492✔
296	x = Vector{keytype(pmf)}(undef, n)	492✔
297	psum = 0.	355✔
298	for (index, p) in pmf	710✔
299	for i in 1:length(x)	2,614✔
300	if psum ≤ u[i] ≤ psum + p	14,916✔
301	x[i] = index	2,151✔
302	end
303	end	27,218✔
304	psum = psum + p	2,614✔
305	end	4,873✔
306	Dict(index => sum(x .== Ref(index)) for index in keys(pmf))	355✔
307	end
308
309	Geometric(p, k) = (1 - p)^k*p	1,836✔
310
311	function update_pmf(estimator::Estimator{<:U})	355✔
312	f(estimator, index) = sqrt(var(estimator, index)/cost(estimator, index))	3,454✔
313	p = interp(f, estimator)	355✔
314	if !any(isnan.(p)) && all(p[2:end] .< 0)	355✔
315	for index in keys(estimator)	154✔
316	if isnan(f(estimator, index)) \|\| isinf(f(estimator, index))	1,672✔
317	set_pmf_key(estimator, index, 2^(p[1]+sum(p[2:end].*index.I)))	628✔
318	else
319	set_pmf_key(estimator, index, f(estimator, index))	522✔
320	end
321	end	1,150✔
322	normalize!(pmf(estimator))	154✔
323	end
324	end
325
326	function normalize!(pmf::Dict{<:Index, Float64})	262✔
327	tot_sum = sum(values(pmf))	524✔
328	for (key, val) in pmf	524✔
329	pmf[key] /= tot_sum	1,906✔
330	end	1,906✔
331	end
332
333	Prob(estimator::Estimator{<:U}, index::Index) = sum(val for (key, val) in pmf(estimator) if key ≥ index)	35,282✔
334
335	Prob(estimator::Estimator{<:U}) = Dict(index => Prob(estimator, index) for index in keys(estimator))	×
336
337	var(estimator::Estimator{<:U, <:MC}, qoi::Integer) = var(accumulator(estimator, qoi))	×
338
339	var(estimator::Estimator{<:U, <:QMC}, qoi::Integer) = mean(var(accumulator(estimator, qoi, shift)) for shift in 1:n_shifts(estimator))	×
340
341	varest(estimator::Estimator{<:U, <:MC}) = aggregate(estimator, map(qoi -> var(accumulator(estimator, qoi)) / length(accumulator(estimator, qoi)), 1:estimator[:nb_of_qoi]))	62,724✔
342
343	varest(estimator::Estimator{<:U, <:QMC}) = aggregate(estimator, map(qoi -> var(mean(accumulator(estimator, qoi, shift)) for shift in 1:n_shifts(estimator)) / n_shifts(estimator), 1:estimator[:nb_of_qoi]))	42,840✔
344
345	n_shifts(estimator::Estimator{<:U, <:QMC}) = size(accumulator(estimator), 2)	79,560✔
346
347	bias(estimator::Estimator{<:U}) = 0.0	1,440✔
348
349	converged(estimator::Estimator{<:U}, ϵ::Real) = mse(estimator) ≤ ϵ^2	×
350
351	mean(estimator::Estimator{<:U}) = aggregate(estimator, map(qoi -> mean(accumulator(estimator, qoi)), 1:estimator[:nb_of_qoi]))	15,120✔

PieterjanRobbe / MultilevelEstimators.jl / 22882410591

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous