From c42aa3d6ed42bf05376a95eba169855d9c6c6e1a Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Mon, 24 Jul 2023 22:37:20 -0400 Subject: [PATCH 1/2] improve docstring --- README.md | 6 ++- docs/src/index.md | 5 +- src/configuration.jl | 82 +++++++++++++++++++++++++++++++-- src/distribution/variable.jl | 6 +-- src/main.jl | 89 +++++++++++++++++++----------------- src/mcmc/montecarlo.jl | 8 ++-- src/statistics.jl | 12 ++--- src/vegas/montecarlo.jl | 4 +- src/vegas_mc/montecarlo.jl | 4 +- 9 files changed, 148 insertions(+), 68 deletions(-) diff --git a/README.md b/README.md index 4ef079f..61c7df7 100644 --- a/README.md +++ b/README.md @@ -24,7 +24,7 @@ julia> f(x, c) = log(x[1]) / sqrt(x[1]) # Define your integrand julia> integrate(f, var = Continuous(0, 1), neval=1e5) # Perform the MC integration for 1e5 steps Integral 1 = -3.99689518016736 ± 0.001364833686666744 (reduced chi2 = 0.695) ``` -In this example, we define an integrand function `f(x, c)` where `x` represents the random variables in the integral and `c` is a `Configuration` object parameter that can hold extra parameters that might be necessary for more complex integrand functions. The variable `x` is determined by the var parameter in `integrate()`. +In this example, we define an integrand function `f(x, c)` where `x` represents the random variables in the integral and `c` is a [`Configuration`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/) object parameter that can hold extra parameters that might be necessary for more complex integrand functions. The variable `x` is determined by the var parameter in `integrate()`. Learn more details from the [documentation](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/). `MCIntegration.jl` also supports Discrete variables. For instance, let's estimate $\pi$ through the Taylor series for $\pi/4 = 1 - 1/3 + 1/5 -1/7 + 1/9 - ...$: ```julia @@ -86,10 +86,12 @@ Here's a brief overview of the three solvers: Given its robustness and efficiency, the default solver in this package is the `:vegasmc`. To choose a specific solver, use the `solver` parameter in the `integrate` function, like `solver=:vegas`. -Please note that the calling convention for the user-defined integrand for `:mcmc` is slightly different from that of `:vegas` and `:vegasmc`. Please refer to the separate detailed note on this. +Please note that the calling convention for the user-defined integrand for `:mcmc` is slightly different from that of `:vegas` and `:vegasmc`. Please refer to the separate [detailed note](https://numericaleft.github.io/MCIntegration.jl/dev/#Algorithm) on this. Packed variables can enhance the efficiency of :vegasmc and :mcmc solvers by reducing the auto-correlation time of the Markov chain, leading to a more effective sampling process. +Learn more from documentation: [Vegas](https://numericaleft.github.io/MCIntegration.jl/dev/lib/vegas/), [VegasMC](https://numericaleft.github.io/MCIntegration.jl/dev/lib/vegasmc/) and [MCMC](https://numericaleft.github.io/MCIntegration.jl/dev/lib/mcmc/) algorithms. + ## Parallelization diff --git a/docs/src/index.md b/docs/src/index.md index dca7169..fff3057 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -196,10 +196,7 @@ julia> plot!(plt, grid, res.mean[2], yerror = res.stdev[2], label="sphere") This package provides three solvers. -- Vegas algorithm (`:vegas`): A Monte Carlo algorithm that uses importance sampling as a variance-reduction technique. Vegas iteratively builds up a piecewise constant weight function, represented -on a rectangular grid. Each iteration consists of a sampling step followed by a refinement -of the grid. The exact details of the algorithm can be found in **_G.P. Lepage, J. Comp. Phys. 27 (1978) 192, 3_** and -**_G.P. Lepage, Report CLNS-80/447, Cornell Univ., Ithaca, N.Y., 1980_**. +- Vegas algorithm (`:vegas`): A Monte Carlo algorithm that uses importance sampling as a variance-reduction technique. Vegas iteratively builds up a piecewise constant weight function, represented on a rectangular grid. Each iteration consists of a sampling step followed by a refinement of the grid. The exact details of the algorithm can be found in **_G.P. Lepage, J. Comp. Phys. 27 (1978) 192, 3_** and **_G.P. Lepage, Report CLNS-80/447, Cornell Univ., Ithaca, N.Y., 1980_**. - Vegas algorithm based on Markov-chain Monte Carlo (`:vegasmc`): A markov-chain Monte Carlo algorithm that uses the Vegas variance-reduction technique. It is as accurate as the vanilla Vegas algorithm, meanwhile tends to be more robust. For complicated high-dimensional integral, the vanilla Vegas algorithm can fail to learn the piecewise constant weight function. This algorithm uses Metropolis–Hastings algorithm to sample the integrand and improves the weight function learning. diff --git a/src/configuration.jl b/src/configuration.jl index 10dfd67..a372dfe 100644 --- a/src/configuration.jl +++ b/src/configuration.jl @@ -36,6 +36,43 @@ 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. """ + +""" + 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 up to now. +- `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. +""" mutable struct Configuration{NI,V,P,O,T} ########### static parameters ################### seed::Int # seed to initialize random numebr generator, also serves as the unique pid of the configuration @@ -91,15 +128,52 @@ By default, dof=[ones(length(var)), ], which means that there is only one integr 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] +- `reweight`: reweight factors for each integrands. If not set, then all factors will be initialized with 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 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. + 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 integrals 1, 2, ..., N and the last entry is for the normalization integral. Only the first integral is connected to the normalization integral. + Only highly correlated integrands 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 """ + +""" + 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 +```julia +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)] +) +``` +""" function Configuration(; var::V=(Continuous(0.0, 1.0),), dof::Union{Int,AbstractVector,AbstractMatrix}=(is_variable(typeof(var)) ? 1 : [ones(Int, length(var)),]), diff --git a/src/distribution/variable.jl b/src/distribution/variable.jl index d3dbb5d..49629e4 100644 --- a/src/distribution/variable.jl +++ b/src/distribution/variable.jl @@ -162,7 +162,7 @@ end """ function Continuous(bounds::AbstractVector{Union{AbstractVector,Tuple}}, size=MaxOrder; offset=0, alpha=2.0, adapt=true, ninc=1000, grid=collect(LinRange(lower, upper, ninc))) -Create a set of continous variable pools sampling from the set [lower, upper) with a distribution generated by a Vegas map, and pack into a ``CompositeVar``. The distribution is trained after each iteraction if `adapt = true`. +Create a set of continous variable pools sampling from the set [lower, upper) with a distribution generated by a Vegas map, and pack into a `CompositeVar`. The distribution is trained after each iteraction if `adapt = true`. # Arguments: - `bounds` : tuple of (lower, upper) for each continuous variable @@ -286,7 +286,7 @@ end """ function Discrete(lower::Int, upper::Int; distribution=nothing, alpha=2.0, adapt=true) -Create a pool of integer variables sampling from the closed set [lower, lower+1, ..., upper] with the distribution `Discrete.distribution``. +Create a pool of integer variables sampling from the closed set [lower, lower+1, ..., upper] with the distribution `Discrete.distribution`. The distribution is trained after each iteraction if `adapt = true`. # Arguments: @@ -330,7 +330,7 @@ end """ function Continuous(bounds::AbstractVector{Union{AbstractVector,Tuple}}, size=MaxOrder; offset=0, alpha=2.0, adapt=true, ninc=1000, grid=collect(LinRange(lower, upper, ninc))) -Create a set of continous variable pools sampling from the set [lower, upper) with a distribution generated by a Vegas map, and pack into a ``CompositeVar``. The distribution is trained after each iteraction if `adapt = true`. +Create a set of continous variable pools sampling from the set [lower, upper) with a distribution generated by a Vegas map, and pack into a `CompositeVar`. The distribution is trained after each iteraction if `adapt = true`. # Arguments: - `bounds` : tuple of (lower, upper) for each continuous variable diff --git a/src/main.jl b/src/main.jl index 256f5bf..cb38c87 100644 --- a/src/main.jl +++ b/src/main.jl @@ -5,58 +5,65 @@ neval=1e4, niter=10, block=16, - verbose=-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, 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 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. - If `inplace = false`, then the signature of the integrand is `integrand(var, config)`, where `var` is a vector of random variables, and `config` is the [`Configuration`](@ref) struct. It should return one or more weights, corresponding to the value of each component of the integrand for the given `var`. - If `inplace = true``, then the signature of the integrand is `integrand(var, weights, config)`, where the additional argument `weights` is the value of the integrand components for the given `var`. - 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`](@ref) 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. -- `verbose`: < -1 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 gamma <=1, where gamma = 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`](@ref), [`VegasMC.montecarlo`](@ref) and [`MCMC.montecarlo`](@ref) 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. -- `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. -- `parallel`: :thread will use Threads.@threads to run different blocks in parallel. Default is :nothread. -- `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`](@ref) docs for more details. +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 ```julia-repl -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 (chi2/dof = 0.645) +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 (chi2/dof = 2.54) +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 integrate(integrand::Function; diff --git a/src/mcmc/montecarlo.jl b/src/mcmc/montecarlo.jl index d0d83cd..5301326 100644 --- a/src/mcmc/montecarlo.jl +++ b/src/mcmc/montecarlo.jl @@ -53,10 +53,10 @@ The third parameter passes the MC `Configuration` struct to the integrand, so th - `measure` : User-defined function with the following signature: ```julia -function measure(idx, var, obs, weight, config) +function measure(idx, var, obs, relative_weight, config) # accumulates the weight into the observable # For example, - # obs[idx] = weight # integral idx + # obs[idx] = relative_weight # integral idx # ... end ``` @@ -80,8 +80,8 @@ The last argument passes the MC `Configuration` struct to the integrand, so that # Examples The following command calls the MC Vegas solver, ```julia-repl -julia> integrate((idx, x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = 2, verbose=-1, solver=:mcmc) -Integral 1 = 0.6757665376867902 ± 0.008655534861083898 (chi2/dof = 0.681) +julia> integrate((idx, x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-1, solver=:mcmc) +Integral 1 = 0.6757665376867902 ± 0.008655534861083898 (reduced chi2 = 0.681) ``` """ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval, diff --git a/src/statistics.jl b/src/statistics.jl index e0e6125..9c78ba4 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -136,16 +136,16 @@ function Base.show(io::IO, ::MIME"text/plain", result::Result) end """ - function report(result::Result, ignore=result.ignore; pick::Union{Function,AbstractVector}=obs -> first(obs), name=nothing, verbose=0) + 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: 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. +- `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. """ function report(result::Result, ignore=result.ignore; pick::Union{Function,AbstractVector}=obs -> first(obs), name=nothing, verbose=0, io::IO=Base.stdout) if isnothing(name) == false @@ -188,7 +188,7 @@ end 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. +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), ...] diff --git a/src/vegas/montecarlo.jl b/src/vegas/montecarlo.jl index 68a61ad..38a8b8e 100644 --- a/src/vegas/montecarlo.jl +++ b/src/vegas/montecarlo.jl @@ -73,8 +73,8 @@ The last argument passes the MC `Configuration` struct to the integrand, so that # Examples The following command calls the Vegas solver, ```julia-repl -julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = 2, verbose=-1, solver=:vegas) -Integral 1 = 0.667203631824444 ± 0.0005046485925614018 (chi2/dof = 1.46) +julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-1, solver=:vegas) +Integral 1 = 0.667203631824444 ± 0.0005046485925614018 (reduced chi2 = 1.46) ``` """ function montecarlo(config::Configuration{Ni,V,P,O,T}, integrand::Function, neval, diff --git a/src/vegas_mc/montecarlo.jl b/src/vegas_mc/montecarlo.jl index 82739e8..7d0404a 100644 --- a/src/vegas_mc/montecarlo.jl +++ b/src/vegas_mc/montecarlo.jl @@ -56,8 +56,8 @@ The last argument passes the MC `Configuration` struct to the integrand, so that # Examples The following command calls the MC Vegas solver, ```julia-repl -julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = 2, verbose=-1, solver=:vegasmc) -Integral 1 = 0.6640840471808533 ± 0.000916060916265263 (chi2/dof = 0.945) +julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-1, solver=:vegasmc) +Integral 1 = 0.6640840471808533 ± 0.000916060916265263 (reduced chi2 = 0.945) ``` """ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval, From aeed69700920d8cc26dde2f53f116c616410e2b1 Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Mon, 24 Jul 2023 23:13:43 -0400 Subject: [PATCH 2/2] improve docstring --- docs/make.jl | 1 + docs/src/lib/utility.md | 5 ++++ src/mcmc/montecarlo.jl | 45 +++++++++++---------------------- src/vegas/montecarlo.jl | 32 +++++++++--------------- src/vegas_mc/montecarlo.jl | 51 +++++++++++++++++++++++++++++++++++++- 5 files changed, 83 insertions(+), 51 deletions(-) create mode 100644 docs/src/lib/utility.md diff --git a/docs/make.jl b/docs/make.jl index f138941..5bc011c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -24,6 +24,7 @@ makedocs(; "lib/vegas.md", "lib/mcmc.md", "lib/distribution.md", + "lib/utility.md", # map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) ] diff --git a/docs/src/lib/utility.md b/docs/src/lib/utility.md new file mode 100644 index 0000000..718c76a --- /dev/null +++ b/docs/src/lib/utility.md @@ -0,0 +1,5 @@ +# Utility + +```@autodocs +Modules = [MCIntegration.MCUtility] +``` \ No newline at end of file diff --git a/src/mcmc/montecarlo.jl b/src/mcmc/montecarlo.jl index 5301326..46540b2 100644 --- a/src/mcmc/montecarlo.jl +++ b/src/mcmc/montecarlo.jl @@ -10,48 +10,38 @@ end verbose=0, timer=[], debug=false; measurefreq::Int=1, measure::Union{Nothing,Function}=nothing, idx::Int=1) where {N,V,P,O,T} -This algorithm calculate high-dimensional integrals with a Markov-chain Monte Carlo. -For multiple integrands invoves multiple variables, it finds the best distribution -ansatz to fit them all together. In additional to the original integral, it also -introduces a normalization integral with integrand ~ 1. +This algorithm computes high-dimensional integrals using a Markov-chain Monte Carlo (MCMC) method. It is particularly well-suited for cases involving multiple integrals over several variables. -Assume we want to calculate the integral ``f_1(x)`` and ``f_2(x, y)``, where x, y are two different types of variables. -The algorithm will try to learn a distribution ``\\rho_x(x)`` and ``\\rho_y(y)`` so that ``f_1(x)/\\rho_x(x)`` and ``f_2(x, y)/\\rho_x(x)/\\rho_y(y)`` -are as flat as possible. +The MCMC algorithm learns an optimal distribution, or 'ansatz', to best fit all integrands under consideration. Additionally, it introduces a normalization integral (with an integrand ~ 1) alongside the original integrals. -The algorithm then samples the variables x and y with a joint distribution using the Metropolis-Hastings algorithm, +Assume we have two integrals to compute, ``f_1(x)`` and ``f_2(x, y)``, where ``x`` and ``y`` are variables of different types. The algorithm aims to learn the distributions ``\\rho_x(x)`` and ``\\rho_y(y)``, such that the quantities ``f_1(x)/\\rho_x(x)`` and ``f_2(x, y)/\\rho_x(x)/\\rho_y(y)`` are as flat as possible. + +Using the Metropolis-Hastings algorithm, the algorithm samples variables ``x`` and ``y`` based on the joint distribution: ```math p(x, y) = r_0 \\cdot \\rho_x(x) \\cdot \\rho_y(y) + r_1 \\cdot |f_1(x)| \\cdot \\rho_y(y) + r_2 \\cdot |f_2(x, y)| ``` -where ``r_i`` are certain reweighting factor to make sure all terms contribute same weights. -One then estimate the integrals by averaging the observables ``f_1(x)\\rho_y(y)/p(x, y)`` and ``f_2(x, y)/p(x, y)``. +where ``r_i`` are reweighting factors ensuring equal contribution from all terms. The integrals are then estimated by averaging the observables ``f_1(x)\\rho_y(y)/p(x, y)`` and ``f_2(x, y)/p(x, y)``. -This algorithm reduces to the vanilla Vegas algorithm by setting ``r_0 = 1`` and ``r_{i>0} = 0``. +Setting ``r_0 = 1`` and ``r_{i>0} = 0`` reduces this algorithm to the classic Vegas algorithm. -The key difference between this algorithmm and the algorithm `:vegasmc` is the way that the joint distribution ``p(x, y)`` is sampled. -In this algorithm, one use Metropolis-Hastings algorithm to sample each term in ``p(x, y)`` as well as the variables ``(x, y)``, so that the MC configuration space consists of -``(idx, x, y)``, where ``idx`` is the index of the user-defined and the normalization integrals. On the other hand, the `:vegasmc` algorithm -only uses Metropolis-Hastings algorithm to sample the configuration space ``(x, y)``. For a given set of x and y, all terms in ``p(x, y)`` are -explicitly calculated on the fly. If one can afford calculating all the integrands on the fly, then `:vegasmc` should be more efficient than this algorithm. +The key difference between this MCMC method and the :vegasmc solver lies in how the joint distribution ``p(x, y)`` is sampled. This MCMC solver uses the Metropolis-Hastings algorithm to sample each term in ``p(x, y)`` as well as the variables ``(x, y)``. The MC configuration space thus consists of `(idx, x, y)`, where `idx` represents the index of the user-defined and normalization integrals. In contrast, the `:vegasmc` algorithm only samples the `(x, y)` space, explicitly calculating all terms in ``p(x, y)`` on-the-fly for a given set of ``x`` and ``y``. -NOTE: If there are more than one integrals, only one of them are sampled and measured at each Markov-chain Monte Carlo step! +Note: When multiple integrals are involved, only one of them is sampled and measured at each Markov-chain Monte Carlo step! -For low-dimensional integrations, this algorithm is much less efficient than the :vegasmc or :vegas solvers. For high-dimension integrations, however, -this algorithm becomes as efficent and robust as the :vegasmc solver, and is more efficient and robust than the :vegas solver. +While the MCMC method may be less efficient than the `:vegasmc or `:vegas` solvers for low-dimensional integrations, it exhibits superior efficiency and robustness when faced with a large number of integrals, a scenario where the `:vegasmc` and `:vegas` solvers tend to struggle. # Arguments -- `integrand` : User-defined function with the following signature: +- `integrand`: A user-defined function evaluating the integrand. The function should be `integrand(idx, var, config)`. Here, `idx` is the index of the integrand component to be evaluated, `var` are the random variables and `weights` is an output array to store the calculated weights. The last parameter passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. + +- `measure`: An optional user-defined function to accumulate the integrand weights into the observable. The function signature should be `measure(idx, var, obs, relative_weight, config)`. Here, `idx` is the integrand index, `obs` is a vector of observable values for each component of the integrand and `relative_weight` is the weight calculated from the `idx`-th integrand multiplied by the probability of the corresponding variables. + +The following are the snippets of the `integrand` and `measure` functions: ```julia function integrand(idx, var, config) # calculate your integrand values # return integrand of the index idx end ``` -The first argument `idx` is index of the integral being sampled. -The second parameter `var` is either a Variable struct if there is only one type of variable, or a tuple of Varibles if there are more than one types of variables. -The third parameter passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. - -- `measure` : User-defined function with the following signature: ```julia function measure(idx, var, obs, relative_weight, config) # accumulates the weight into the observable @@ -60,11 +50,6 @@ function measure(idx, var, obs, relative_weight, config) # ... end ``` -The first argument `idx` is index of the integral being sampled. -The second argument `var` is either a Variable struct if there is only one type of variable, or a tuple of Varibles if there are more than one types of variables. -The third argument passes the user-defined observable to the function, it should be a vector with the length same as the integral number. -The fourth argument is the integrand weights to be accumulated to the observable, it is a vector with the length same as the integral number. -The last argument passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. # Remark: diff --git a/src/vegas/montecarlo.jl b/src/vegas/montecarlo.jl index 38a8b8e..9bfeea2 100644 --- a/src/vegas/montecarlo.jl +++ b/src/vegas/montecarlo.jl @@ -24,38 +24,34 @@ end # end """ - function montecarlo(config::Configuration{Ni,V,P,O,T}, integrand::Function, neval, verbose=0, timer=[], debug=false; measure::Union{Nothing,Function}=nothing, measurefreq::Int=1) where {Ni,V,P,O,T} -This algorithm implements the classic Vegas algorithm. - -The main idea of the algorithm can be found in this [link](https://vegas.readthedocs.io/en/latest/background.html#importance-sampling). +This function implements the Vegas algorithm, a Monte Carlo method specifically designed for multi-dimensional integration. The underlying principle and methodology of the algorithm can be explored further in the [Vegas documentation](https://vegas.readthedocs.io/en/latest/background.html#importance-sampling). -The algorithm uses a simple important sampling scheme. -Consider an one-dimensional integral with the integrand ``f(x)``, the algorithm will try to learn an -optimized distribution ``\\rho(x)>0`` which mimic the integrand as good as possible (a.k.a, the Vegas map trick, see [`Dist.Continuous`](@ref)) for more detail. +# Overview +The Vegas algorithm employs an importance sampling scheme. For a one-dimensional integral with the integrand ``f(x)``, the algorithm constructs an optimized distribution ``\\rho(x)`` that approximates the integrand as closely as possible (a strategy known as the Vegas map trick; refer to [`Dist.Continuous`](@ref) for more details). -One then sample the variable ``x`` with the distribution ``\\rho(x)``, and estimate the integral by averging the estimator ``f(x)/\\rho(x)``. +The variable ``x`` is then sampled using the distribution ``\\rho(x)``, and the integral is estimated by averaging the estimator ``f(x)/\\rho(x)``. -NOTE: If there are more than one integrals, then all integrals are sampled and measured at each Monte Carlo step. +# Note +- If there are multiple integrals, all of them are sampled and measured at each Monte Carlo step. +- This algorithm is particularly efficient for low-dimensional integrations but might be less efficient and robust than the Markov-chain Monte Carlo solvers for high-dimensional integrations. -This algorithm is very efficient for low-dimensional integrations, but can be less -efficient and less robust than the Markov-chain Monte Carlo solvers for high-dimensional integrations. # Arguments -- `integrand` : User-defined function with the following signature: +- `integrand`: A user-defined function evaluating the integrand. 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. The last parameter passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. + +- `measure`: An optional user-defined function to accumulate the integrand weights into the observable. 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. + +The following are the snippets of the `integrand` and `measure` functions: ```julia function integrand(var, config) # calculate your integrand values # return integrand1, integrand2, ... end ``` -The first parameter `var` is either a Variable struct if there is only one type of variable, or a tuple of Varibles if there are more than one types of variables. -The second parameter passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. - -- `measure` : User-defined function with the following signature: ```julia function measure(var, obs, weights, config) # accumulates the weight into the observable @@ -65,10 +61,6 @@ function measure(var, obs, weights, config) # ... end ``` -The first argument `var` is either a Variable struct if there is only one type of variable, or a tuple of Varibles if there are more than one types of variables. -The second argument passes the user-defined observable to the function, it should be a vector with the length same as the integral number. -The third argument is the integrand weights to be accumulated to the observable, it is a vector with the length same as the integral number. -The last argument passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. # Examples The following command calls the Vegas solver, diff --git a/src/vegas_mc/montecarlo.jl b/src/vegas_mc/montecarlo.jl index 7d0404a..6bad860 100644 --- a/src/vegas_mc/montecarlo.jl +++ b/src/vegas_mc/montecarlo.jl @@ -1,5 +1,4 @@ """ - function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval, verbose=0, debug=false; measurefreq::Int=1, measure::Union{Nothing,Function}=nothing) where {N,V,P,O,T} @@ -60,6 +59,56 @@ julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2, Integral 1 = 0.6640840471808533 ± 0.000916060916265263 (reduced chi2 = 0.945) ``` """ + +""" + function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval, + verbose=0, debug=false; + measurefreq::Int=1, measure::Union{Nothing,Function}=nothing) where {N,V,P,O,T} + +This function applies a Markov-chain Monte Carlo (MCMC) technique combined with the Vegas algorithm to compute integrals. In addition to calculating the original integrals, it also introduces a normalization integral with an integrand ~ 1, which enhances the efficiency and robustness of high-dimensional integration tasks. + +# Overview +Given multiple integrands involving multiple variables, the algorithm finds the best distribution ansatz that fits all integrands together. For instance, consider we want to calculate two integrals: ``f_1(x)`` and ``f_2(x, y)``, where ``x`` and ``y`` are two different types of variables. The algorithm learns distributions ``\\rho_x(x)`` and ``\\rho_y(y)`` such that ``f_1(x)/\\rho_x(x)`` and ``f_2(x, y)/\\rho_x(x)/\\rho_y(y)`` are as flat as possible. + +Then, it samples variables ``x`` and ``y`` using the Metropolis-Hastings algorithm with a joint distribution `p(x, y)`, +```math +p(x, y) = r_0 \\cdot \\rho_x(x) \\cdot \\rho_y(y) + r_1 \\cdot |f_1(x)| \\cdot \\rho_y(y) + r_2 \\cdot |f_2(x, y)| +``` +where ``r_i`` are certain reweighting factor to make sure all terms contribute same weights. + +One can then estimate the integrals by averaging the observables ``f_1(x)\\rho_y(y)/p(x, y)`` and ``f_2(x, y)/p(x, y)``. + +The algorithm defaults to the standard Vegas algorithm if ``r_0 = 1`` and ``r_{i>0} = 0``. + +# Arguments +- `integrand`: A user-defined function evaluating the integrand. 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. The last parameter passes the MC `Configuration` struct to the integrand, so that user has access to userdata, etc. + +- `measure`: An optional user-defined function to accumulate the integrand weights into the observable. 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. + +The following are the snippets of the `integrand` and `measure` functions: +```julia +function integrand(var, config) + # calculate your integrand values + # return integrand1, integrand2, ... +end +``` +```julia +function measure(var, obs, weights, config) + # accumulates the weight into the observable + # For example, + # obs[1] = weights[1] # integral 1 + # obs[2] = weights[2] # integral 2 + # ... +end +``` + +# Examples +The following command calls the VegasMC solver, +```julia-repl +julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-1, solver=:vegasmc) +Integral 1 = 0.6640840471808533 ± 0.000916060916265263 (reduced chi2 = 0.945) +``` +""" function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval, verbose=0, timer=[], debug=false; measure::Union{Nothing,Function}=nothing, measurefreq::Int=1, inplace::Bool=false