diff --git a/.gitignore b/.gitignore index 6e24bde..074f0af 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ Manifest.toml /docs/build/ +.*vscode/ \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 8c9481b..9de690a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -4,24 +4,24 @@ using Documenter DocMeta.setdocmeta!(MetropolisAlgorithm, :DocTestSetup, :(using MetropolisAlgorithm); recursive = true) makedocs(; - modules = [MetropolisAlgorithm], - authors = "Shuhei Ohno", - sitename = "MetropolisAlgorithm.jl", - format = Documenter.HTML(; - canonical = "https://JuliaFewBody.github.io/MetropolisAlgorithm.jl", - edit_link = "main", - assets = String[ - "./assets/logo.ico", + modules = [MetropolisAlgorithm], + authors = "Shuhei Ohno", + sitename = "MetropolisAlgorithm.jl", + format = Documenter.HTML(; + canonical = "https://JuliaFewBody.github.io/MetropolisAlgorithm.jl", + edit_link = "main", + assets = String[ + "./assets/logo.ico", + ], + ), + pages = [ + "Home" => "index.md", + "Examples" => "examples.md", + "API reference" => "API.md", ], - ), - pages = [ - "Home" => "index.md", - "Examples" => "examples.md", - "API reference" => "API.md", - ], ) deploydocs(; - repo = "github.com/JuliaFewBody/MetropolisAlgorithm.jl", - devbranch = "main", + repo = "github.com/JuliaFewBody/MetropolisAlgorithm.jl", + devbranch = "main", ) diff --git a/docs/src/examples.md b/docs/src/examples.md index ba20bd6..6586a39 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -34,7 +34,7 @@ X = [r[1] for r in R] hist(X) ``` -This is the trajectory of a walker at each step. The histogram above shows the number of these dots in the interval of bins. +This is the trajectory of a walker at each step. The histogram above shows the number of these points in each bin. ```@example eg1 # plotting trajectory diff --git a/docs/src/index.md b/docs/src/index.md index ee0e043..f23a086 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -35,7 +35,7 @@ r₀ = [0.0, 0.0] # initial position R = metropolis(p, r₀) # 1-walker x 10000-steps ``` -In the variational Monte Carlo method (VMC), sampling is performed simultaneously with multiple walkers (rather than just one walker). Here is an example of 10000-walkers x 5-steps, 2-dimensional Metropolis-walk using `metropolis!()`. This function overwrites its second argument without memory allocation, where the first argument is the (normalized or unnormalized) distribution function, and the second argument is the vector of the initial value vectors. Use the For statement to repeat as many times as you like. Remove the first several steps by yourself to ensure equilibrium. +In the variational Monte Carlo method (VMC), sampling is performed simultaneously with multiple walkers (rather than just one walker). Here is an example of 10000-walkers x 5-steps, 2-dimensional Metropolis-walk using `metropolis!()`. This function overwrites its second argument without memory allocation, where the first argument is the (normalized or unnormalized) distribution function, and the second argument is the vector of the initial value vectors. Use the For statement to repeat as many times as you like. Discard the first several steps to ensure equilibrium. ```@example using MetropolisAlgorithm diff --git a/src/MetropolisAlgorithm.jl b/src/MetropolisAlgorithm.jl index 350256f..7beea02 100644 --- a/src/MetropolisAlgorithm.jl +++ b/src/MetropolisAlgorithm.jl @@ -9,136 +9,137 @@ Random.seed!(123) # one-step x many-walkers function metropolis!(f::Function, R::Vector{<:Vector}; type = typeof(first(first(R))), d::Real = one(type)) - # initialize - half = type(1//2) - n_dim = length(first(R)) - # Metropolis-walk - for i ∈ keys(R) - # shift - Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ - r_old = R[i] - f_old = f(r_old) - r_new = r_old .+ Δr - f_new = f(r_new) - # accept - p = min(1, f_new / f_old) - if rand() < p - R[i] = r_new + # initialize + half = type(1 // 2) + n_dim = length(first(R)) + # Metropolis-walk + for i in keys(R) + # shift + Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ + r_old = R[i] + f_old = f(r_old) + r_new = r_old .+ Δr + f_new = f(r_new) + # accept + p = min(1, f_new / f_old) + if rand() < p + R[i] = r_new + end end - end - return + return end # many-steps x one-walker function metropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}; type = typeof(first(r_ini)), d::Real = one(type)) - # initialize - half = type(1//2) - n_dim = length(r_ini) - n_steps = length(R) - r_old = r_ini - f_old = f(r_ini) - # Metropolis-walk - for i ∈ 2:n_steps - # shift - Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ - r_new = r_old + Δr - f_new = f(r_new) - # accept - p = min(1, f_new / f_old) - if rand() < p - r_old = r_new - f_old = f_new + # initialize + half = type(1 // 2) + n_dim = length(r_ini) + n_steps = length(R) + r_old = r_ini + f_old = f(r_ini) + # Metropolis-walk + for i in 2:n_steps + # shift + Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ + r_new = r_old + Δr + f_new = f(r_new) + # accept + p = min(1, f_new / f_old) + if rand() < p + r_old = r_new + f_old = f_new + end + # save + R[i] = r_old end - # save - R[i] = r_old - end - return + return end # many-steps x one-walker with memory allocation function metropolis(f::Function, r_ini::Vector{<:Real}; n_steps::Int = 10^4, type = typeof(first(r_ini)), d::Real = one(type)) - # memory allocation - R = fill(typeof(r_ini)(undef, size(r_ini)), n_steps) - R[begin] = r_ini - # Metropolis sampling - metropolis!(f, R, r_ini; d = d) - # return - return R + # memory allocation + R = fill(typeof(r_ini)(undef, size(r_ini)), n_steps) + R[begin] = r_ini + # Metropolis sampling + metropolis!(f, R, r_ini; d = d) + # return + return R end struct bin - min::Vector{<:Real} - max::Vector{<:Real} - width::Vector{<:Real} - number::Vector{<:Int} - center::Vector{Vector{<:Real}} - corner::Vector{Vector{<:Real}} - counter::Array - function bin(A::Vector{<:Vector}; number = fill(10, length(first(A)))) # A = [[0,2], [2,2], [2,4]], number = [3,3] - min = [minimum(a[i] for a ∈ A) for i ∈ keys(first(A))] # [0,2] - max = [maximum(a[i] for a ∈ A) for i ∈ keys(first(A))] # [2,4] - width = [(max[n] - min[n]) / (number[n]-1) for n ∈ keys(number)] # [1,1] - center = [[((i)*max[n] + (number[n]-1-i)*min[n]) / (number[n]-1) for i ∈ 0:(number[n] - 1)] for n ∈ keys(number)] # [[0,1,2], [2,3,4]] - corner = [[((i-1//2)*max[n] + (number[n]-1-i+1//2)*min[n]) / (number[n]-1) for i ∈ 0:number[n]] for n ∈ keys(number)] # [[-0.5,0.5,1.5,2.5], [1.5,2.5,3.5,4.5]] - counter = zeros(Int64, number...) - for k ∈ Iterators.product([1:n for n ∈ number]...) - a₋ = [corner[i][k[i]] for i ∈ keys(k)] - a₊ = [corner[i][k[i] + 1] for i ∈ keys(k)] - c = count(a -> prod(a₋ .≤ a .< a₊), A) - counter[k...] = c + min::Vector{<:Real} + max::Vector{<:Real} + width::Vector{<:Real} + number::Vector{<:Int} + center::Vector{Vector{<:Real}} + corner::Vector{Vector{<:Real}} + counter::Array + function bin(A::Vector{<:Vector}; number = fill(10, length(first(A)))) # A = [[0,2], [2,2], [2,4]], number = [3,3] + min = [minimum(a[i] for a in A) for i in keys(first(A))] # [0,2] + max = [maximum(a[i] for a in A) for i in keys(first(A))] # [2,4] + width = [(max[n] - min[n]) / (number[n] - 1) for n in keys(number)] # [1,1] + center = [[((i) * max[n] + (number[n] - 1 - i) * min[n]) / (number[n] - 1) for i in 0:(number[n] - 1)] for n in keys(number)] # [[0,1,2], [2,3,4]] + corner = [[((i - 1 // 2) * max[n] + (number[n] - 1 - i + 1 // 2) * min[n]) / (number[n] - 1) for i in 0:number[n]] for n in keys(number)] # [[-0.5,0.5,1.5,2.5], [1.5,2.5,3.5,4.5]] + counter = zeros(Int64, number...) + for k in Iterators.product([1:n for n in number]...) + a₋ = [corner[i][k[i]] for i in keys(k)] + a₊ = [corner[i][k[i] + 1] for i in keys(k)] + c = count(a -> prod(a₋ .≤ a .< a₊), A) + counter[k...] = c + end + return new(min, max, width, number, center, corner, counter) end - new(min, max, width, number, center, corner, counter) - end end function Base.count(center::Vector, width::Vector, A::Vector{<:Vector}) - a₋ = center .- width / 2 - a₊ = center .+ width / 2 - return count(a -> prod(a₋ .≤ a .< a₊), A) + a₋ = center .- width / 2 + a₊ = center .+ width / 2 + return count(a -> prod(a₋ .≤ a .< a₊), A) end function pdf(center::Vector, width::Vector, A::Vector{<:Vector}) - return count(center, width, A) / length(A) / prod(width) + return count(center, width, A) / length(A) / prod(width) end # docstrings @doc raw""" - metropolis!(f::Function, R::Vector{<:Vector}; type=typeof(first(first(R))), d::Real=one(type)) + metropolis!(f::Function, R::Vector{<:Vector}; type=typeof(first(first(R))), d::Real=one(type)) -This function calculates **one-step of many-walkers** and overwrites the second argument `R`. Each child vector in `R` is a point of the walker (not a trajectory). +This function performs **one step for many walkers** and overwrites the second argument `R`. +Each element of `R` is a point (not a trajectory). # Arguments - `f::Function`: Distribution function. It does not need to be normalized. -- `R::Vector{<:Vector}`: Vector of vectors (points). Each child vector is a point of the walker. -- `type::Type=typeof(first(first(R)))`: Type of trajectory points. e.g., Float32, Float64, etc.. -- `d::Real=one(type)`: Maximum step size. +- `R::Vector{<:Vector}`: Vector of vectors (points). Each element is a point of a walker. +- `type::Type=typeof(first(first(R)))`: Type of trajectory points, e.g. Float32 or Float64. +- `d::Real=one(type)`: Maximum step size. """ metropolis!(f::Function, R::Vector{<:Vector}) @doc raw""" - metropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}; type=typeof(first(r_ini)), d::Real=one(type)) + metropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}; type=typeof(first(r_ini)), d::Real=one(type)) -This function calculates **many-steps of one-walker** and overwrites the second argument `R`. Each child vector in `R` is a point of the trajectory. +This function performs **many steps for one walker** and overwrites the second argument `R`. # Arguments - `f::Function`: Distribution function. It does not need to be normalized. -- `r_ini::Vector{<:Real}`: Initial value vector. Even in the one-dimensional case, the initial value must be defined as a vector. Each child vector (point) has the same size as `r_ini`. -- `R::Vector{<:Vector}`: Vector of vectors (points). Each child vector is a point of the walker. The first element of `R` is same as `r_ini`. -- `n_steps::Int=10^5`: Number of steps. It is same as the length of the output parent vector. -- `type::Type=typeof(first(r_ini))`: Type of trajectory points. e.g., Float32, Float64, etc.. +- `r_ini::Vector{<:Real}`: Initial value vector. Even in the one-dimensional case the initial value must be a vector. Each element (point) has the same size as `r_ini`. +- `R::Vector{<:Vector}`: Vector of vectors (points). Each element is a point of the trajectory. The first element of `R` is the same as `r_ini`. +- `n_steps::Int=10^4`: Number of steps; this is the length of the output vector `R` and matches the default used in `metropolis`. +- `type::Type=typeof(first(r_ini))`: Type of trajectory points, e.g. Float32 or Float64. - `d::Real=one(type)`: Maximum step size. Default value is 1. """ metropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}) @doc raw""" - metropolis(f::Function, r_ini::Vector{<:Real}; n_steps::Int=10^5, type=typeof(first(r_ini)), d::Real=one(type)) + metropolis(f::Function, r_ini::Vector{<:Real}; n_steps::Int=10^4, type=typeof(first(r_ini)), d::Real=one(type)) -This function calculates **many-steps of one-walker** using `metropolis!(f, R, r_ini)` and returns the trajectory `R` as a vector of vectors (points) with memory allocation. +This function performs **many steps for one walker** using `metropolis!(f, R, r_ini)` and returns the trajectory `R` as a vector of vectors (points) with memory allocation. """ metropolis(f::Function, r_ini::Vector{<:Real}) @doc raw""" - bin(A::Vector{<:Vector}; number = fill(10,length(first(A)))) + bin(A::Vector{<:Vector}; number = fill(10,length(first(A)))) -This function creates a data for multidimensional histogram for testing. +This function creates data for a multidimensional histogram (for testing). # Examples ```julia @@ -199,7 +200,6 @@ julia> exp(0) / sqrt(2*π) julia> pdf([0.0, 0.0], [0.2, 0.2], [randn(2) for i in 1:1000000]) 0.15389999999999998 -s julia> exp(0) / sqrt(2*π)^2 0.15915494309189537