Manual

Index

IndexSet

MultilevelEstimators.FTType
FT(d::Integer)
FT(δ::Real...)
FT(δ::AbstractVector{<:Real})
FT(δ::NTuple{N, <:Real})

Return a full tensor index set in d dimenions with optional weights δ. Default weights are all 1's.

Examples

julia> FT(2)
FT{2}

julia> FT([1, 1/2, 1/2])
FT{3}

julia> print(FT(2), 4)
  ◼ ◼ ◼ ◼ ◼
  ◼ ◼ ◼ ◼ ◼
  ◼ ◼ ◼ ◼ ◼
  ◼ ◼ ◼ ◼ ◼
  ◼ ◼ ◼ ◼ ◼

See also: SL, ML, TD, HC, ZC, AD, U

source
MultilevelEstimators.TDType
TD(d::Integer)
TD(δ::Real...)
TD(δ::AbstractVector{<:Real})
TD(δ::NTuple{N, <:Real})

Return a total degree index set in d dimenions with optional weights δ. Default weights are all 1's.

Examples

julia> TD(2)
TD{2}

julia> TD([1, 1/2, 1/2])
TD{3}

julia> print(TD(2), 4)
  ◼
  ◼ ◼
  ◼ ◼ ◼
  ◼ ◼ ◼ ◼
  ◼ ◼ ◼ ◼ ◼

See also: SL, ML, FT, HC, ZC, AD, U

source
MultilevelEstimators.HCType
HC(d::Integer)
HC(δ::Real...)
HC(δ::AbstractVector{<:Real})
HC(δ::NTuple{N, <:Real})

Return a hyperbolic cross index set in d dimenions with optional weights δ. Default weights are all 1's.

Examples

julia> HC(2)
HC{2}

julia> HC([1, 1/2, 1/2])
HC{3}

julia> print(HC(2), 4)
  ◼
  ◼
  ◼
  ◼ ◼
  ◼ ◼ ◼ ◼ ◼

See also: SL, ML, FT, TD, ZC, AD, U

source
MultilevelEstimators.ZCType
ZC(d::Integer)
ZC(δ::Real...)
ZC(δ::AbstractVector{<:Real})
ZC(δ::NTuple{N, <:Real})

Return a Zaremba cross index set in d dimenions with optional weights δ. Default weights are all 1's.

Examples

julia> ZC(2)
ZC{2}

julia> ZC([1, 1/2, 1/2])
ZC{3}

julia> print(ZC(2), 4)
  ◼ ◼ 
  ◼ ◼ 
  ◼ ◼ ◼ 
  ◼ ◼ ◼ ◼ ◼ 
  ◼ ◼ ◼ ◼ ◼ 

See also: SL, ML, FT, TD, HC, AD, U

source
MultilevelEstimators.get_index_setFunction
get_index_set(idxset::AbstractIndexSet, sz::Integer)

Return an iterator over all indices in the index set idxset for a given size parameter sz.

Examples

julia> collect(get_index_set(TD(2), 2))
6-element Vector{CartesianIndex{2}}:
 CartesianIndex(0, 0)
 CartesianIndex(1, 0)
 CartesianIndex(2, 0)
 CartesianIndex(0, 1)
 CartesianIndex(1, 1)
 CartesianIndex(0, 2)
source
Base.printFunction
print(idxset::AbstractIndexSet{2}, sz::Integer)

Print the indices in the index set idxset for a given size parameter sz. Only implemented for two-dimensional index sets.

Examples

julia> print(TD(2), 3)
  ◼
  ◼ ◼
  ◼ ◼ ◼
  ◼ ◼ ◼ ◼
source

SampleMethod

Distribution

MultilevelEstimators.transformFunction
transform(D::AbstractDistribution, x::Real)

Apply a transformation to the uniformly distributed number x such that it is distributed according to D.

The transform is implemented by applying the inverse CDF of the distribution D.

Note

This is currently limited only to distributions with an analytic expression for the inverse CDF.

Examples

julia> x = rand(10)
10-element Vector{Float64}:
 0.0491718221481211
 0.11907881640750706
 0.3932710232252806
 0.024094310524527707
 0.6918572875342215
 0.7675180540873912
 0.08725304891274233
 0.8557176841095734
 0.8025607099234905
 0.661425351684768

julia> broadcast(i -> transform(Normal(), i), x)
10-element Vector{Float64}:
 -1.6529372045666144
 -1.1796042977687324
 -0.27080365779705895
 -1.975701269628406
  0.5011217721531482
  0.7306975836397528
 -1.3578663461442697
  1.0612756865840536
  0.8508033865054415
  0.41635630569809523
source

LatticeRule32

MultilevelEstimators.LatticeRule32Type
LatticeRule32(z, s, n)
LatticeRule32(z, s)
LatticeRule32(z)

Returns a rank-1 lattice rule in s dimensions with generating vector z and at most n points.

When no maximum number of points n is provided, we assume n=2^32. When no number of dimensions s is provided, we assume s=length(z).

Info

Technically, we return an extensible lattice sequence where the k-th point is transformed using the radical inverse function. This has the advantage that we can add points to the set without changing the already computed points.

More generating vectors can be found online here or here.

Examples

julia> lattice_rule = LatticeRule32([0x00000001,0x00000022], 2, 0x00000037) # Fibonacci lattice rule
LatticeRule32{2}

julia> get_point(lattice_rule, 2)
2-element Vector{Float64}:
 0.75
 0.5

julia> get_point(lattice_rule, 56) # returns random points beyond n
2-element Vector{Float64}:
 0.07336635446929285
 0.34924148955718615

See also: get_point, ShiftedLatticeRule

source
LatticeRule32(str, s, n)
LatticeRule32(str, s)
LatticeRule32(str)

Returns a rank-1 lattice rule in s dimensions with generating vector z read from the file str and at most n points.

Examples

julia> z_file = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "K_3600_32.txt");

julia> lattice_rule = LatticeRule32(z_file, 16)
LatticeRule32{16}

julia> get_point(lattice_rule, 123)
16-element Vector{Float64}:
 0.3828125
 0.2890625
 0.1484375
 0.3515625
 0.6015625
 0.1171875
 0.4296875
 0.4609375
 0.2265625
 0.6484375
 0.9609375
 0.6796875
 0.0546875
 0.2265625
 0.1328125
 0.1640625

See also: get_point, ShiftedLatticeRule

source
MultilevelEstimators.ShiftedLatticeRuleType
ShiftedLatticeRule(lat::AbstractLatticeRule)

Returns a shifted rank-1 lattice rule based on the lattice rule lat.

Examples

julia> z_file = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "K_3600_32.txt");

julia> lattice_rule = LatticeRule32(z_file, 16)
LatticeRule32{16}

julia> shifted_lattice_rule = ShiftedLatticeRule(lattice_rule)
ShiftedLatticeRule{LatticeRule32{16}}

julia> get_point(shifted_lattice_rule, 0)
16-element Vector{Float64}:
 0.0491718221481211
 0.11907881640750706
 0.3932710232252806
 0.024094310524527707
 0.6918572875342215
 0.7675180540873912
 0.08725304891274233
 0.8557176841095734
 0.8403841370820818
 0.8907696748195567
 0.138227024723224
 0.3477368477058109
 0.19852090963358837
 0.008197858843363659
 0.5920407556271403
 0.8010551449181886

See also: LatticeRule32, get_point

source
MultilevelEstimators.get_pointFunction
get_point(lat::AbstractLatticeRule, k::Integer)

Get the k-th point of the lattice rule lat.

julia> lattice_rule = LatticeRule32([0x00000001,0x00000022], 2, 0x00000037) # Fibonacci lattice rule
LatticeRule32{2}

julia> get_point(lattice_rule, 2)
2-element Vector{Float64}:
 0.75
 0.5

See also: LatticeRule32, ShiftedLatticeRule

source

Estimator

MultilevelEstimators.EstimatorType
Estimator(index_set::AbstractIndexSet,
          sample_method::AbstractSampleMethod,
          sample_function::Function,
          distributions::Vector{<:Distribution}; kwargs...)

Create an Estimator with index set index_set and sample method sample_method to compute the expected value of the quantity of interest returned by sample_function, and where distributions is the uncertainty on the input parameters.

Examples

See the example section in the online documentation for examples on how to use this function.

General keywords

Different algorithmic options can be set by providing appropriate keyword arguments. A list of these keyword arguments and their description is provided below.

Note

Not all keyword arguments can be used in combination with all Estimator types.

name – a human readable name for the problem that must be solved. This name will be used when saving a History file with diagnostic information from each run, and also when generating reports with the Reporter-package. Default is UntitledEstimator.

folder is where the History files with diagnostic information from every run will be saved. Default is the current directory.

save – if set to true, MultilevelEstimators saves diagnostic information about the simulation in a separate file in the .jld2 format. See JLD2 for more information about this file format. By default, this keyword is set to true.

verbose – if set to true, MultilevelEstimators will print useful information about the run to STDOUT. Default is true.

nb_of_warm_up_samples is the number of warm-up samples to take at each level (or index) to obtain initial variance estimates. In combination with do_regression = true (default), this is only true for all levels < 2 (or indices with all coordinate directions < 2). On all other levels (or indices), this is the initial maximum number of regressed samples. When using MC-sampling, this value must be larger than or equal to 2, in order to be able to compute the variance. The default value is 20. When using QMC-sampling, this is the number of samples per random shift, and since we only consider the variance with repect to all random shifts, this value can be as small as 1 (default). Increase this number for more accurate initial variance estimation.

nb_of_qoi is the number of quantities of interest returned by the sample_function. In the default case (aggregation = "max"), the simulation will be performed by controlling the root mean square error of the quantity of interest with the largest variance. The expected value of all other quantities of interest is guaranteed to be estimated more accurate or as accurate as the one with the largest variance. The default value is 1. The quantitiy of interest can optionally be provided using the key qoi_with_max_var, to ensure that the same quantity is used in every iteration. Other modes of aggregation are "sum" and "norm" (see aggregation).

aggregation specifies the way in which the variance from multiple quantities of interest must be aggregated to determine the optimal number of samples. Options are max, sum or norm. The default value is max.

continuate – if set to true, each call to run with a specified tolerance tol results in a simulation for a sequence of larger tolerances $\text{tol}_i = p^{(n-1-i)} \text{tol}$ to get improved estimates for the variance and the bias. The default value is true.

continuation_mul_factor is the value $p$ in the above formula. The default value is 1.2.

nb_of_tols is the value $n$ in the above formula. The default value is 10.

save_samples – if set to true, all samples of a run will be saved together with the diagnostic information. This is usefull for plotting histograms etc.

cost_model allows the user to provide a function that returns the cost of a sample on each index. For example, a geometric cost model could be provided with

γ = 1.5
Estimator(...; cost_model = level -> 2^(γ * level[1]))

When no cost model is provided (default), the actual run time on each level (or index) is used as a proxy for the cost model.

Tip

Beware of precompilation! Using actual run times for the cost model might be very unreliable, especially on the coarsest level (or index), when the Estimator is run for the first time. In real production code, it is advised that the run-function is called once beforehand with a larger tolerance than the target accuracuy, to eliminate precompilation time.

run(estimator, 10 * tol) # precompilation run
run(estimator, tol)

nb_of_workers allows the user to specify the number of workers (processors) to be used for parallel sampling in a multicore environment. By default, this is equal to the value returned by nworkers(). nb_of_workers can also be given as function, specifying the number of workers on each level (or index). This is particularly usefull for expensive sample functions that require non-trivial load balancing.

nb_of_uncertainties is a function indicating the number of uncertainties to be used on each level (or each index). The default value is equal to the length of distributions, the vector of Distributions specifying the input uncertainty of the parameters. This is useful for so-called level-dependent estimators, where a different number of parameters is activated on each level (or index). The length of distributions must be larger than or equal to the maximum value returned by this function.

max_index_set_param is the maximum index set size parameter that can be used in a simulation. Typically, a simulation will only be possible up to a finite index set size parameter, due to, for instance, memory constraints or prohibitively lare computational cost. For a multilevel simulation with ML this is the maximum number of levels. For a multi-index simulation with FT, TD, HC or ZC, this is the maximum index set size parameter. For simulations with AD or U, this is the maximum index set size to be used in the max_search_space (see below). Default is 10 times the number of dimensions of the index set. When this value is exceeded, the simulation will stop with a bias (and, consequently, a root mean square error) larger than the requested accuracy.

min_index_set_param is the minimum index set size parameter that should be used in a simulation. Useful for providing a minimum number of levels that should be used. This parameter should be at least 2 to allow for regression and estimation of the rates alpha and `beta. Default value is 2.

Keywords valid for index set types SL, ML, FT, TD, HC and ZC

do_regression – if set to true, the variance on the finer levels (or indices) will be guessed from the aready available variances on the coarser levels (or indices). These finer levels (or indices) must have a value larger than or equal to 2 in at least one of the coordinate directions. This is very powerfull when dealing with computationally demanding sample functions. Default: true.

do_mse_splitting – if set to true, MultilevelEstimators will use a non-trivial splitting of the mean square error of the estimator: $\text{MSE} = \theta \text{var} + (1 - \theta) \text{bias}^2$. The value of the splitting parameter is very important to obtain an efficient estimator. Default is true.

min_splitting is the minimum value for the mean square error splitting parameter. Default: 0.5.

max_splitting is the maximum value for the mean square error splitting parameter. Default: 0.99.

Keywords valid for index set types AD and U

max_search_space determines, together with max_index_set_param, the collection of all indices that may be considered for further refinement. This is necessary for AD- and U-type index sets, that do not define such a set, by defintion. Default: TD.

Keywords valid for AD index sets

penalization is a factor added to the computation of the profit indicators in the adaptive algorithm. We observed that in some cases, the algorithm stretches its search space too much around the coordinate axes, especially when actual computation times are used (see cost_model). The penalization parameter 0<p<1 penalizes these search directions, in favor of more classically-shaped index sets, such as TD.

\[P_\ell = \frac{E_\ell}{(\sqrt{V_\ell W_\ell})^p}\]

accept_rate is the acceptance rate in the accept-reject method for chosing new indices from the active set for further refinement. A lower value means a more global search stretegy. Usefull when the adaptive method finds itself stuck along a single coordinate direction. Default value is 1 (no accept-reject).

Keywords valid for QMC sample methods

nb_of_shifts indicates the number of randomly-shifted rank-1 lattice rules used in the QMC method. The higher this value, the more accurate the variance estimation will be. This keyword can also be a function that returns the number of shifts on each level (or index). Typically, the number of shifts on coarser levels will be smaller than the number of shifts on finer levels, since, on these coarse levels, many QMC samples will be available. Default is 10 random shifts.

Note

We recommend a value for nb_of_shifts between 10 and 30.

sample_mul_factor is the multplication factor in the QMC "doubling algorithm". Typically, and as can be justified theoretically, this value is equal to 2. We find that it is sometimes more beneficial to use a smaller multiplication factor, since direclty doubling the number of samples is a quite dramatic event. Default: 1.2.

point_generator is a keyword that allows the user to specify a custom QMC point set of type AbstractLatticeRule. Default is a 3600-dimensional lattice rule from this website.

See also: run, History

source
Base.runFunction

run(estimator::Estimator, ε::Real) run(estimator::Estimator, ε::Vector{<:Real})

Run the estimator and compute the expected value of the quantity of interest up to the given tolerance(s) ε. Returns a History-object that contains usefull diagnostics about the simulation.

Examples

An example detailing how to use MultilevelEstimators is provided in the example in the documentation.

See also: Estimator, History

source

History

MultilevelEstimators.HistoryType
History

A History is the result of a call to run and contains useful diagnostic information about the state of the Estimator.

The information is stored as a Vector of Dict's, where each entry in the vector is the result of a single run. See below for a list of keys added to these Dicts. Access to the key corresponding to the last entry in the History (that is, the one with the highest accuracy) is provided as history[key].

See also the package Reporter for automatic generation of reports based on a History object.

Keys

  • type::AbstractIndexSet: the index set type of the Estimator.
  • ndims::Integer: the number of dimensions of the Estimator.
  • name::String: the name of this Estimator.
  • folder::String: the directory where the results are stored.
  • elapsed::Real: the elapsed time (in seconds).
  • tol::Real: the request accuracy.
  • current_index_set::Vector{::Index}: the levels (or indices) currently in use by the Estimator.
  • index_set::Vector{::Index}: all levels (or indices) that are added to the Estimator.
  • mse::Real: the measured mean square error of the Estimator.
  • rmse::Real: the measured root mean square error of the Estimator.
  • mean::Real: the estimated expected value of the quantity of interest, up to an absolute error of tol.
  • var::Real: the measured variance of the quantity of interest.
  • varest::Real: the measured variance of the Estimator for the expected value of the quantity of interest.
  • bias::Real: the measured bias of the Estimator.
  • E::Dict{::Index, ::Real}: a Dict with the measured expected value of the quantity of interest on each level (or index).
  • V::Dict{::Index, ::Real}: a Dict with the measured variance of the quantity of interest on each level (or index).
  • dE::Dict{::Index, ::Real}: a Dict with the measured expected value of the difference of the quantity of interest on each level (or index).
  • dV::Dict{::Index, ::Real}: a Dict with the measured variance of the difference of the quantity of interest on each level (or index).
  • T::Dict{::Index, ::Real}: a Dict with the measured computation time for a sample of the quantity of interest on each level (or index).
  • W::Dict{::Index, ::Real}: a Dict with the computational cost for a sample of the quantity of interest on each level (or index). Only available when a cost model was provided!
  • α::Tuple: the estimated rates of decay of the expected value of the difference of the quantity of interest.
  • β::Tuple: the estimated rates of decay of the variance of the difference of the quantity of interest.
  • γ::Tuple: the estimated rates of increase of the cost of the difference of the quantity of interest.
  • nb_of_samples::Array: an Array with, on each index, the number of samples taken on that index.
  • logbook: a logbook with information about the adaptive algorithm Only available for Estimators with index sets of type AD
  • samples::Array: an Array with on each index, the collection of samples of the quantity of interest taken on that index. Only available when the optional key save_samples = true was used
  • samples_diff::Array: an Array with on each index, the collection of samples of the difference of the quantity of interest taken on that index. Only available when the optional key save_samples = true was used

See also: Estimator, run

source