Monte Carlo Integrator

mutable struct Configuration

Struct that contains everything needed for MC.

Static parameters

• seed: seed to initialize random numebr generator, also serves as the unique pid of the configuration
• rng: a MersenneTwister random number generator, seeded by seed
• userdata: user-defined parameter
• var: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.

integrand properties

• neighbor::Vector{Tuple{Int, Int}} : vector of tuples that defines the neighboring integrands. Two neighboring integrands are directly connected in the Markov chain. e.g., [(1, 2), (2, 3)] means the integrand 1 and 2 are neighbor, and 2 and 3 are neighbor. The neighbor vector defines a undirected graph showing how the integrands are connected. Please make sure all integrands are connected. By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to [(N+1, 1), (1, 2), (2, 4), ..., (N-1, N)], where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram. Only highly correlated integrands are not highly correlated should be defined as neighbors. Otherwise, most of the updates between the neighboring integrands will be rejected and wasted.
• dof::Vector{Vector{Int}}: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2.
• observable: observables that is required to calculate the integrands, will be used in the measure function call. It is either an array of any type with the common operations like +-*/^ defined.
• reweight: reweight factors for each integrands. The reweight factor of the normalization diagram is assumed to be 1. Note that you don't need to explicitly add the normalization diagram.
• visited: how many times this integrand is visited by the Markov chain.

current MC state

• step: the number of MC updates performed up to now
• norm: the index of the normalization diagram. norm is larger than the index of any user-defined integrands
• normalization: the accumulated normalization factor. Physical observable = Configuration.observable/Configuration.normalization.
• propose/accept: array to store the proposed and accepted updates for each integrands and variables. Their shapes are (number of updates X integrand number X max(integrand number, variable number). The last index will waste some memory, but the dimension is small anyway.

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(type, length(dof)),
    reweight::Vector{Float64}=ones(length(dof) + 1),
    seed::Int=rand(Random.RandomDevice(), 1:1000000),
    neighbor::Union{Vector{Vector{Int}},Vector{Tuple{Int,Int}},Nothing}=nothing,
    userdata=nothing,
    kwargs...
)

Create a Configuration struct

Arguments

• var: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.

By default, var = (Continuous(0.0, 1.0),), which is a single continuous variable.

• dof::Vector{Vector{Int}}: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2.

By default, dof=[ones(length(var)), ], which means that there is only one integrand, and each variable has one degree of freedom.

• obs: observables that is required to calculate the integrands, will be used in the measure function call.

It is either an array of any type with the common operations like +-*/^ defined. By default, it will be set to 0.0 if there is only one integrand (e.g., length(dof)==1); otherwise, it will be set to zeros(length(dof)).

• para: user-defined parameter, set to nothing if not needed
• reweight: reweight factors for each integrands. If not set, then all factors will be initialized with one.
• seed: seed to initialize random numebr generator, also serves as the unique pid of the configuration. If it is nothing, then use RandomDevice() to generate a random seed in [1, 1000_1000]
• neighbor::Vector{Tuple{Int, Int}} : vector of tuples that defines the neighboring integrands. Two neighboring integrands are directly connected in the Markov chain. e.g., [(1, 2), (2, 3)] means the integrand 1 and 2 are neighbor, and 2 and 3 are neighbor. The neighbor vector defines a undirected graph showing how the integrands are connected. Please make sure all integrands are connected. By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to [(N+1, 1), (1, 2), (2, 4), ..., (N-1, N)], where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram. Only highly correlated integrands are not highly correlated should be defined as neighbors. Otherwise, most of the updates between the neighboring integrands will be rejected and wasted.
• userdata: User data you want to pass to the integrand and the measurement

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
• chi2: chi-square per dof of the MC integration
• neval: number of evaluations of the integrand
• ignore: ignore iterations untill ignore
• dof: degrees of freedom of the MC integration (number of iterations - 1)
• 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, 
    print=-1, 
    gamma=1.0, 
    adapt=true,
    debug=false, 
    reweight_goal::Union{Vector{Float64},Nothing}=nothing, 
    ignore::Int=adapt ? 1 : 0,
    measure::Union{Nothing,Function}=nothing,
    measurefreq::Int=1,
    kwargs...
)

Calculate the integrals, collect statistics, and return a Result struct that contains the estimations and errors.

Remarks

• User may run the MC in parallel using MPI. Simply run mpiexec -n N julia userscript.jl where N is the number of workers. In this mode, only the root process returns meaningful results. All other workers return nothing, nothing. User is responsible to handle the returning results properly. If you have multiple number of mpi version, you can use "mpiexecjl" in your "~/.julia/package/MPI/###/bin" to make sure the version is correct. See https://juliaparallel.github.io/MPI.jl/stable/configuration/ for more detail.
• In the MC, a normalization diagram is introduced to normalize the MC estimates of the integrands. More information can be found in the link: https://kunyuan.github.io/QuantumStatistics.jl/dev/man/important_sampling/#Important-Sampling. User don't need to explicitly specify this normalization diagram.Internally, normalization diagram will be added to each table that is related to the integrands.

Arguments

• integrand:Function call to evaluate the integrand. It should accept an argument of the type Configuration, and return a weight. Internally, MC only samples the absolute value of the weight. Therefore, it is also important to define Main.abs for the weight if its type is user-defined.
• solver : :vegas, :vegasmc, or :mcmc. See Readme for more details.
• config: Configuration object to perform the MC integration. If nothing, it attempts to create a new one with Configuration(; kwargs...).
• neval: Number of evaluations of the integrand per iteration.
• niter: Number of iterations. The reweight factor and the variables will be self-adapted after each iteration.
• block: Number of blocks. Each block will be evaluated by about neval/block times. Each block is assumed to be statistically independent, and will be used to estimate the error. In MPI mode, the blocks are distributed among the workers. If the numebr of workers N is larger than block, then block will be set to be N.
• print: -2 to not print anything; -1 to print minimal information; 0 to print the iteration history in the end; >0 to print MC configuration for every print seconds and print the iteration history in the end.
• gamma: Learning rate of the reweight factor after each iteraction. Note that alpha <=1, where alpha = 0 means no reweighting.
• adapt: Whether to adapt the grid and the reweight factor.
• debug: Whether to print debug information (type instability, float overflow etc.)
• reweight_goal: The expected distribution of visited times for each integrand after reweighting . If not set, then all factors will be initialized with one. Only useful for the :mcmc solver.
• ignore: ignore the iteration until the ignore round. By default, the first iteration is igonred if adapt=true, and non is ignored if adapt=false.
• measure: measurement function, See Vegas.montecarlo, VegasMC.montecarlo and MCMC.montecarlo for more details.
• measurefreq: how often perform the measurement for ever measurefreq MC steps. If a measurement is expansive, you may want to make the measurement less frequent.
• kwargs: Keyword arguments. The supported keywords include,
  • measure and measurefreq: measurement function and how frequent it is called.
  • If config is nothing, you may need to provide arguments for the Configuration constructor, check Configuration docs for more details.

Examples

julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = 2, print=-1)
Integral 1 = 0.6668625385256122 ± 0.0009960142738129768   (chi2/dof = 1.07)

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

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: the ignore the first # iteractions.
• 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.