You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
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.
38
+
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.
Copy file name to clipboardExpand all lines: src/MetropolisAlgorithm.jl
+86-86Lines changed: 86 additions & 86 deletions
Original file line number
Diff line number
Diff line change
@@ -9,136 +9,137 @@ Random.seed!(123)
9
9
10
10
# one-step x many-walkers
11
11
functionmetropolis!(f::Function, R::Vector{<:Vector}; type =typeof(first(first(R))), d::Real=one(type))
12
-
# initialize
13
-
half =type(1//2)
14
-
n_dim =length(first(R))
15
-
# Metropolis-walk
16
-
for i ∈keys(R)
17
-
# shift
18
-
Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ
19
-
r_old = R[i]
20
-
f_old =f(r_old)
21
-
r_new = r_old .+ Δr
22
-
f_new =f(r_new)
23
-
# accept
24
-
p =min(1, f_new / f_old)
25
-
ifrand() < p
26
-
R[i] = r_new
12
+
# initialize
13
+
half =type(1//2)
14
+
n_dim =length(first(R))
15
+
# Metropolis-walk
16
+
for i inkeys(R)
17
+
# shift
18
+
Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ
19
+
r_old = R[i]
20
+
f_old =f(r_old)
21
+
r_new = r_old .+ Δr
22
+
f_new =f(r_new)
23
+
# accept
24
+
p =min(1, f_new / f_old)
25
+
ifrand() < p
26
+
R[i] = r_new
27
+
end
27
28
end
28
-
end
29
-
return
29
+
return
30
30
end
31
31
32
32
# many-steps x one-walker
33
33
functionmetropolis!(f::Function, R::Vector{<:Vector}, r_ini::Vector{<:Real}; type =typeof(first(r_ini)), d::Real=one(type))
34
-
# initialize
35
-
half =type(1//2)
36
-
n_dim =length(r_ini)
37
-
n_steps =length(R)
38
-
r_old = r_ini
39
-
f_old =f(r_ini)
40
-
# Metropolis-walk
41
-
for i ∈2:n_steps
42
-
# shift
43
-
Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ
44
-
r_new = r_old + Δr
45
-
f_new =f(r_new)
46
-
# accept
47
-
p =min(1, f_new / f_old)
48
-
ifrand() < p
49
-
r_old = r_new
50
-
f_old = f_new
34
+
# initialize
35
+
half =type(1//2)
36
+
n_dim =length(r_ini)
37
+
n_steps =length(R)
38
+
r_old = r_ini
39
+
f_old =f(r_ini)
40
+
# Metropolis-walk
41
+
for i in2:n_steps
42
+
# shift
43
+
Δr = d * (rand(type, n_dim) .- half) # [-d/2,d/2)ⁿ
44
+
r_new = r_old + Δr
45
+
f_new =f(r_new)
46
+
# accept
47
+
p =min(1, f_new / f_old)
48
+
ifrand() < p
49
+
r_old = r_new
50
+
f_old = f_new
51
+
end
52
+
# save
53
+
R[i] = r_old
51
54
end
52
-
# save
53
-
R[i] = r_old
54
-
end
55
-
return
55
+
return
56
56
end
57
57
58
58
# many-steps x one-walker with memory allocation
59
59
functionmetropolis(f::Function, r_ini::Vector{<:Real}; n_steps::Int=10^4, type =typeof(first(r_ini)), d::Real=one(type))
60
-
# memory allocation
61
-
R =fill(typeof(r_ini)(undef, size(r_ini)), n_steps)
62
-
R[begin] = r_ini
63
-
# Metropolis sampling
64
-
metropolis!(f, R, r_ini; d = d)
65
-
# return
66
-
return R
60
+
# memory allocation
61
+
R =fill(typeof(r_ini)(undef, size(r_ini)), n_steps)
62
+
R[begin] = r_ini
63
+
# Metropolis sampling
64
+
metropolis!(f, R, r_ini; d = d)
65
+
# return
66
+
return R
67
67
end
68
68
69
69
struct bin
70
-
min::Vector{<:Real}
71
-
max::Vector{<:Real}
72
-
width::Vector{<:Real}
73
-
number::Vector{<:Int}
74
-
center::Vector{Vector{<:Real}}
75
-
corner::Vector{Vector{<:Real}}
76
-
counter::Array
77
-
functionbin(A::Vector{<:Vector}; number =fill(10, length(first(A)))) # A = [[0,2], [2,2], [2,4]], number = [3,3]
78
-
min = [minimum(a[i] for a ∈ A) for i ∈keys(first(A))] # [0,2]
79
-
max = [maximum(a[i] for a ∈ A) for i ∈keys(first(A))] # [2,4]
80
-
width = [(max[n] - min[n]) / (number[n]-1) for n ∈keys(number)] # [1,1]
81
-
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]]
82
-
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]]
83
-
counter =zeros(Int64, number...)
84
-
for k ∈ Iterators.product([1:n for n ∈ number]...)
85
-
a₋ = [corner[i][k[i]] for i ∈keys(k)]
86
-
a₊ = [corner[i][k[i] +1] for i ∈keys(k)]
87
-
c =count(a ->prod(a₋ .≤ a .< a₊), A)
88
-
counter[k...] = c
70
+
min::Vector{<:Real}
71
+
max::Vector{<:Real}
72
+
width::Vector{<:Real}
73
+
number::Vector{<:Int}
74
+
center::Vector{Vector{<:Real}}
75
+
corner::Vector{Vector{<:Real}}
76
+
counter::Array
77
+
functionbin(A::Vector{<:Vector}; number =fill(10, length(first(A)))) # A = [[0,2], [2,2], [2,4]], number = [3,3]
78
+
min = [minimum(a[i] for a in A) for i inkeys(first(A))] # [0,2]
79
+
max = [maximum(a[i] for a in A) for i inkeys(first(A))] # [2,4]
80
+
width = [(max[n] - min[n]) / (number[n] -1) for n inkeys(number)] # [1,1]
81
+
center = [[((i) * max[n] + (number[n] -1- i) * min[n]) / (number[n] -1) for i in0:(number[n] -1)] for n inkeys(number)] # [[0,1,2], [2,3,4]]
82
+
corner = [[((i -1//2) * max[n] + (number[n] -1- i +1//2) * min[n]) / (number[n] -1) for i in0:number[n]] for n inkeys(number)] # [[-0.5,0.5,1.5,2.5], [1.5,2.5,3.5,4.5]]
83
+
counter =zeros(Int64, number...)
84
+
for k in Iterators.product([1:n for n in number]...)
85
+
a₋ = [corner[i][k[i]] for i inkeys(k)]
86
+
a₊ = [corner[i][k[i] +1] for i inkeys(k)]
87
+
c =count(a ->prod(a₋ .≤ a .< a₊), A)
88
+
counter[k...] = c
89
+
end
90
+
returnnew(min, max, width, number, center, corner, counter)
89
91
end
90
-
new(min, max, width, number, center, corner, counter)
91
-
end
92
92
end
93
93
94
94
function Base.count(center::Vector, width::Vector, A::Vector{<:Vector})
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).
109
+
This function performs **one step for many walkers** and overwrites the second argument `R`.
110
+
Each element of `R` is a point (not a trajectory).
110
111
111
112
# Arguments
112
113
- `f::Function`: Distribution function. It does not need to be normalized.
113
-
- `R::Vector{<:Vector}`: Vector of vectors (points). Each child vector is a point of the walker.
114
-
- `type::Type=typeof(first(first(R)))`: Type of trajectory points. e.g., Float32, Float64, etc..
115
-
- `d::Real=one(type)`: Maximum step size.
114
+
- `R::Vector{<:Vector}`: Vector of vectors (points). Each element is a point of a walker.
115
+
- `type::Type=typeof(first(first(R)))`: Type of trajectory points, e.g. Float32 or Float64.
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.
122
+
This function performs **manysteps for onewalker** and overwrites the second argument `R`.
122
123
123
124
# Arguments
124
125
- `f::Function`: Distribution function. It does not need to be normalized.
125
-
- `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`.
126
-
- `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`.
127
-
- `n_steps::Int=10^5`: Number of steps. It is same as the length of the output parent vector.
128
-
- `type::Type=typeof(first(r_ini))`: Type of trajectory points. e.g., Float32, Float64, etc..
126
+
- `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`.
127
+
- `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`.
128
+
- `n_steps::Int=10^4`: Number of steps; this is the length of the output vector `R` and matches the default used in `metropolis`.
129
+
- `type::Type=typeof(first(r_ini))`: Type of trajectory points, e.g. Float32 or Float64.
129
130
- `d::Real=one(type)`: Maximum step size. Default value is 1.
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.
136
+
This function performs **manysteps for onewalker** using `metropolis!(f, R, r_ini)` and returns the trajectory `R` as a vector of vectors (points) with memory allocation.
136
137
"""metropolis(f::Function, r_ini::Vector{<:Real})
137
138
138
139
@docraw"""
139
-
bin(A::Vector{<:Vector}; number = fill(10,length(first(A))))
140
+
bin(A::Vector{<:Vector}; number = fill(10,length(first(A))))
140
141
141
-
This function creates a data for multidimensional histogram for testing.
142
+
This function creates data for a multidimensional histogram (for testing).
142
143
143
144
# Examples
144
145
```julia
@@ -199,7 +200,6 @@ julia> exp(0) / sqrt(2*π)
199
200
200
201
julia> pdf([0.0, 0.0], [0.2, 0.2], [randn(2) for i in 1:1000000])
0 commit comments