Main module

mutable struct Configuration{NI,V,P,O,T}

Struct that holds all the necessary parameters and variables for Monte Carlo integration.

Parameters

• NI : Number of integrands
• V : Type of variables
• P : Type of user-defined data
• O : Type of observables
• T : Type of integrand

Static parameters

• seed: Seed to initialize the random number generator, also serves as the unique process ID of the configuration.
• rng: A MersenneTwister random number generator, seeded by seed.
• userdata: User-defined parameter.
• var: Tuple of variables. Each variable should derive from the abstract type Variable (see variable.jl for details). Using a tuple instead of a vector improves performance.

Integrands properties

• neighbor::Vector{Tuple{Int, Int}} : Vector of tuples defining neighboring integrands. Two neighboring integrands are directly connected in the Markov chain. The neighbor vector defines an undirected graph showing how the integrands are connected. Only highly correlated integrands should be defined as neighbors to reduce autocorrelations.
• dof::Vector{Vector{Int}}: Degrees of freedom of each integrand, i.e., the dimensions in which each integrand can vary.
• observable: Observables required to calculate the integrands, will be used in the measure function call.
• reweight: Reweight factors for each integrand. The reweight factor of the normalization integrand (namely, the last element) is assumed to be 1.
• visited: The number of times each integrand is visited by the Markov chain.

Current MC state

• step: The number of Monte Carlo updates performed thus far.
• norm: The index of the normalization integrand. norm is larger than the index of any user-defined integrands.
• normalization: The accumulated normalization factor.
• propose/accept: Arrays to store the proposed and accepted updates for each integrand and variable.

function Configuration(;
    var::Union{Variable,AbstractVector,Tuple}=(Continuous(0.0, 1.0),),
    dof::Union{Int,AbstractVector,AbstractMatrix}=[ones(Int, length(var))],
    type=Float64,  # type of the integrand
    obs::AbstractVector=zeros(Float64, length(dof)),
    reweight::Vector{Float64}=ones(length(dof) + 1),
    seed::Int=rand(Random.RandomDevice(), 1:1000000),
    userdata=nothing,
    neighbor::Union{Vector{Vector{Int}},Vector{Tuple{Int,Int}},Nothing}=nothing,
    kwargs...
)

Create a Configuration struct for MC integration.

Arguments

• var: Either a single Variable, a CompositeVar, or a tuple consisting of Variable and/or CompositeVar instances. Tuples are used to improve performance by ensuring type stability. By default, var is set to a tuple containing a single continuous variable, (Continuous(0.0, 1.0),).
• dof: Degrees of freedom for each integrand, as a vector of integers. For example, [[0, 1], [2, 3]] means the first integrand has zero instances of var#1 and one of var#2; while the second integrand has two instances of var#1 and 3 of var#2. Defaults to [ones(length(var))], i.e., one degree of freedom for each variable.
• type: Type of the integrand, Float64 by default.
• obs: Vector of observables needed to calculate the integrands, which will be used in the measure function call.
• reweight: Vector of reweight factors for each integrand. By default, all factors are initialized to one. Internally, a reweight factor of 1 will be appended to the end of the reweight vector, which is for the normalization integral.
• seed: Seed for the random number generator. This also serves as the unique identifier of the configuration. If it is nothing, then a random seed between 1 and 1,000,000 is generated.
• userdata: User data to pass to the integrand and the measurement.
• neighbor: Vector of tuples that define neighboring integrands. For example, [(1, 2), (2, 3)] means that the first and second integrands, and the second and third integrands, are neighbors. Neighboring integrands are directly connected in the Markov chain. By default, all integrands are connected in ascending order. Note that the normalization integral is automatically appended at the end of the integrand list and is considered as neighbor with the first user-defined integrand.

Example

cfg = Configuration(
    var = (Continuous(0.0, 1.0), Continuous(-1.0, 1.0)),
    dof = [[1, 1], [2, 0]],
    obs = [0.0, 0.0],
    seed = 1234,
    neighbor = [(1, 2)]
)

struct Result{O,C}

the returned result of the MC integration.

Members

• mean: mean of the MC integration
• stdev: standard deviation of the MC integration samples
• chi2: reduced chi-square of the MC integration samples
• neval: number of evaluations of the integrand
• ignore: ignore iterations untill ignore
• config: configuration of the MC integration from the last iteration
• iterations: list of tuples [(data, error, Configuration), ...] from each iteration

function average(history, idx=1; init=1, max=length(history))

Average the history[1:max]. Return the mean, standard deviation and chi2 of the history.

Arguments

• history: a list of tuples, such as [(data, error, Configuration), ...]
• idx: the index of the integral
• max: the last index of the history to average with
• init : the first index of the history to average with

function integrate(integrand::Function;
    solver::Symbol=:vegas, # :mcmc, :vegas, or :vegasmc
    config::Union{Configuration,Nothing}=nothing,
    neval=1e4, 
    niter=10, 
    block=16, 
    measure::Union{Nothing,Function}=nothing,
    measurefreq::Int=1,
    inplace::Bool=false,
    adapt=true,
    gamma=1.0, 
    reweight_goal::Union{Vector{Float64},Nothing}=nothing, 
    parallel::Symbol=:nothread,
    ignore::Int=adapt ? 1 : 0,
    debug=false, 
    verbose=-1, 
    kwargs...
)

Calculate the integrals, collect statistics, and return a Result struct containing the estimates and errors.

Arguments

• integrand: A user-provided function to compute the integrand values. The function signature differs based on the selected solver and whether computations are done in-place:
  • For solver = :vegas or :vegasmc, the function should be either integrand(var, config) or integrand(var, weights, config) depending on whether inplace is false or true respectively. Here, var are the random variables and weights is an output array to store the calculated weights.
  • For solver = :mcmc, the function should be integrand(idx, var, config), where idx is the index of the integrand component to be evaluated.

Keyword Arguments

• solver: Integration algorithm to use: :vegas, :vegasmc, or :mcmc. Default is :vegas.
• config: Configuration object for the integration. If nothing, a new one is created using Configuration(; kwargs...).
• neval: Number of integrand evaluations per iteration (default: 1e4).
• niter: Number of iterations for the integration process (default: 10).
• block: Number of blocks for statistical independence assumption (default: 16).
• measure: An optional measurement function.
  • For solver = :vegas or :vegasmc, the function signature should be measure(var, obs, relative_weights, config). Here, obs is a vector of observable values for each component of the integrand and relative_weights are the weights calculated from the integrand multiplied by the probability of the corresponding variables.
  • For solver = :mcmc, the signature should be measure(idx, var, obs, relative_weight, config), where obs is the observable vector and relative_weight is the weight calculated from the idx-th integrand multiplied by the probability of the variables.
• measurefreq: How often the measurement function is called (default: 1).
• inplace: Whether to use the inplace version of the integrand. Default is false, which is more convenient for integrand with a few return values but may cause type instability. Only useful for the :vegas and :vegasmc solver.
• adapt: Whether to adapt the grid and the reweight factor (default: true).
• gamma: Learning rate of the reweight factor after each iteration (default: 1.0).
• reweight_goal: The expected distribution of visited times for each integrand after reweighting. Default is nothing.
• parallel: Run different blocks in parallel. Options are :thread and :nothread. Default is :nothread.
• ignore: Ignore the iteration until the ignore round. By default, the first iteration is ignored if adapt=true, and none is ignored if adapt=false.
• verbose: Control the printing level of the iteration history and configuration.
  • <-1:print nothing
  • -1: print minimal information (Default)
  • 0: print iteration history
  • >0: print MC configuration every verbose seconds and print iteration history.
• debug: Whether to print debug information such as type instability or float overflow (default: false).
• kwargs: Other keyword arguments for the Configuration constructor.

Returns

Returns a Result struct containing the estimates and errors of the calculated integrals.

Notes

• In MPI mode, only the root process returns meaningful results. All other workers return nothing. Users should handle the returning results properly.

• The solvers :vegasmc and :vegas automatically append a normalization integral to the end of the integrand vector. When providing reweight_goal, don't forget assign the weight (the last element) for this normalization integral.

Examples

integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-2, solver=:vegas)
Integral 1 = 0.6663652080622751 ± 0.000490978424216832   (reduced chi2 = 0.645)

julia> integrate((x, f, c)-> (f[1] = x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-2, solver=:vegas, inplace=true)
Integral 1 = 0.6672083165915914 ± 0.0004919147870306026   (reduced chi2 = 2.54)

function report(result::Result, ignore=result.ignore; pick::Union{Function,AbstractVector}=obs -> first(obs), name=nothing, verbose=0, io::IO=Base.stdout)

print the summary of the result. It will first print the configuration from the last iteration, then print the weighted average and standard deviation of the picked observable from each iteration.

Arguments

• result: Result object contains the history from each iteration
• ignore: ignore the first # iterations.
• pick: The pick function is used to select one of the observable to be printed. The return value of pick function must be a Number.
• name: name of each picked observable. If name is not given, the index of the pick function will be used.