From cf81504e993766adb952a78980c880540eb97582 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Mon, 8 Sep 2025 21:18:07 +0200 Subject: [PATCH 1/3] updated wording and formatting --- .gitignore | 1 + docs/make.jl | 32 +++---- src/MetropolisAlgorithm.jl | 172 ++++++++++++++++++------------------- 3 files changed, 103 insertions(+), 102 deletions(-) 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/src/MetropolisAlgorithm.jl b/src/MetropolisAlgorithm.jl index 350256f..74b1a5f 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 From 65deb5247b55cdbb9968602f08a79fd57331a610 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Mon, 8 Sep 2025 21:30:40 +0200 Subject: [PATCH 2/3] minor typos in docs and formatting --- docs/src/examples.md | 2 +- docs/src/index.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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 From b1c373c2d3f085fb800b800f4d6d50fcd830573f Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 9 Sep 2025 17:53:05 +0200 Subject: [PATCH 3/3] updated to bold font --- src/MetropolisAlgorithm.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MetropolisAlgorithm.jl b/src/MetropolisAlgorithm.jl index 74b1a5f..7beea02 100644 --- a/src/MetropolisAlgorithm.jl +++ b/src/MetropolisAlgorithm.jl @@ -106,7 +106,7 @@ end @doc raw""" metropolis!(f::Function, R::Vector{<:Vector}; type=typeof(first(first(R))), d::Real=one(type)) -This function performs one step for many walkers and overwrites the second argument `R`. +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 @@ -119,7 +119,7 @@ Each element of `R` is a point (not a trajectory). @doc raw""" metropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}; type=typeof(first(r_ini)), d::Real=one(type)) -This function performs many steps for one walker and overwrites the second argument `R`. +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. @@ -133,7 +133,7 @@ This function performs many steps for one walker and overwrites the second argum @doc raw""" metropolis(f::Function, r_ini::Vector{<:Real}; n_steps::Int=10^4, type=typeof(first(r_ini)), d::Real=one(type)) -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. +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"""