Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
Manifest.toml
/docs/build/
.*vscode/
32 changes: 16 additions & 16 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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",
)
2 changes: 1 addition & 1 deletion docs/src/examples.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
172 changes: 86 additions & 86 deletions src/MetropolisAlgorithm.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
Loading