|
1 | | -adapt_step_size{T<:Hamiltonian}(spl::Sampler{T}, stats::Float64, δ::Float64) = begin |
| 1 | +type WarmUpManager |
| 2 | + iter_n :: Int |
| 3 | + state :: Int |
| 4 | + params :: Dict |
| 5 | +end |
| 6 | + |
| 7 | +getindex(wum::WarmUpManager, param) = wum.params[param] |
| 8 | + |
| 9 | +setindex!(wum::WarmUpManager, value, param) = wum.params[param] = value |
| 10 | + |
| 11 | +init_warm_up_params{T<:Hamiltonian}(vi::VarInfo, spl::Sampler{T}) = begin |
| 12 | + wum = WarmUpManager(1, 1, Dict()) |
| 13 | + |
| 14 | + # Pre-cond |
| 15 | + wum[:θ_num] = 0 |
| 16 | + wum[:θ_mean] = nothing |
| 17 | + D = length(vi[spl]) |
| 18 | + wum[:stds] = ones(D) |
| 19 | + wum[:vars] = ones(D) |
| 20 | + |
| 21 | + # DA |
| 22 | + wum[:ϵ] = nothing |
| 23 | + wum[:μ] = nothing |
| 24 | + wum[:ϵ_bar] = 1.0 |
| 25 | + wum[:H_bar] = 0.0 |
| 26 | + wum[:m] = 0 |
| 27 | + wum[:n_adapt] = spl.alg.n_adapt |
| 28 | + wum[:δ] = spl.alg.delta |
| 29 | + |
| 30 | + spl.info[:wum] = wum |
| 31 | +end |
| 32 | + |
| 33 | +update_da_params(wum::WarmUpManager, ϵ::Float64) = begin |
| 34 | + wum[:ϵ] = [ϵ] |
| 35 | + wum[:μ] = log(10 * ϵ) |
| 36 | +end |
| 37 | + |
| 38 | +adapt_step_size(wum::WarmUpManager, stats::Float64) = begin |
2 | 39 | dprintln(2, "adapting step size ϵ...") |
3 | | - m = spl.info[:m] += 1 |
4 | | - if m <= spl.alg.n_adapt |
5 | | - γ = 0.05; t_0 = 10; κ = 0.75 |
6 | | - μ = spl.info[:μ]; ϵ_bar = spl.info[:ϵ_bar]; H_bar = spl.info[:H_bar] |
| 40 | + m = wum[:m] += 1 |
| 41 | + if m <= wum[:n_adapt] |
| 42 | + γ = 0.05; t_0 = 10; κ = 0.75; δ = wum[:δ] |
| 43 | + μ = wum[:μ]; ϵ_bar = wum[:ϵ_bar]; H_bar = wum[:H_bar] |
7 | 44 |
|
8 | 45 | H_bar = (1 - 1 / (m + t_0)) * H_bar + 1 / (m + t_0) * (δ - stats) |
9 | 46 | ϵ = exp(μ - sqrt(m) / γ * H_bar) |
10 | 47 | dprintln(1, " ϵ = $ϵ, stats = $stats") |
11 | 48 |
|
12 | 49 | ϵ_bar = exp(m^(-κ) * log(ϵ) + (1 - m^(-κ)) * log(ϵ_bar)) |
13 | | - push!(spl.info[:ϵ], ϵ) |
14 | | - spl.info[:ϵ_bar], spl.info[:H_bar] = ϵ_bar, H_bar |
| 50 | + push!(wum[:ϵ], ϵ) |
| 51 | + wum[:ϵ_bar], wum[:H_bar] = ϵ_bar, H_bar |
15 | 52 |
|
16 | | - if m == spl.alg.n_adapt |
| 53 | + if m == wum[:n_adapt] |
17 | 54 | dprintln(0, " Adapted ϵ = $ϵ, $m HMC iterations is used for adaption.") |
18 | 55 | end |
19 | 56 | end |
20 | 57 | end |
21 | 58 |
|
22 | | -init_da_parameters{T<:Hamiltonian}(spl::Sampler{T}, ϵ::Float64) = begin |
23 | | - spl.info[:ϵ] = [ϵ] |
24 | | - spl.info[:μ] = log(10 * ϵ) |
25 | | - spl.info[:ϵ_bar] = 1.0 |
26 | | - spl.info[:H_bar] = 0.0 |
27 | | - spl.info[:m] = 0 |
28 | | -end |
| 59 | +update_pre_cond(wum::WarmUpManager, θ_new) = begin |
| 60 | + |
| 61 | + wum[:θ_num] += 1 # θ_new = x_t |
| 62 | + t = wum[:θ_num] # t |
29 | 63 |
|
30 | | -update_pre_cond{T<:Hamiltonian}(vi::VarInfo, spl::Sampler{T}) = begin |
31 | | - θ_new = realpart(vi[spl]) # x_t |
32 | | - spl.info[:θ_num] += 1 |
33 | | - t = spl.info[:θ_num] # t |
34 | | - θ_mean_old = copy(spl.info[:θ_mean]) # x_bar_t-1 |
35 | | - spl.info[:θ_mean] = (t - 1) / t * spl.info[:θ_mean] + θ_new / t # x_bar_t |
36 | | - θ_mean_new = spl.info[:θ_mean] # x_bar_t |
37 | | - |
38 | | - if t == 2 |
39 | | - first_two = [θ_mean_old'; θ_new'] # θ_mean_old here only contains the first θ |
40 | | - spl.info[:θ_vars] = diag(cov(first_two)) |
41 | | - elseif t <= 1000 |
42 | | - D = length(θ_new) |
43 | | - # D = 2.4^2 |
44 | | - spl.info[:θ_vars] = (t - 1) / t * spl.info[:θ_vars] .+ 100 * eps(Float64) + |
45 | | - (2.4^2 / D) / t * (t * θ_mean_old .* θ_mean_old - (t + 1) * θ_mean_new .* θ_mean_new + θ_new .* θ_new) |
| 64 | + if t == 1 |
| 65 | + wum[:θ_mean] = θ_new |
| 66 | + else |
| 67 | + θ_mean_old = copy(wum[:θ_mean]) # x_bar_t-1 |
| 68 | + wum[:θ_mean] = (t - 1) / t * wum[:θ_mean] + θ_new / t # x_bar_t |
| 69 | + θ_mean_new = wum[:θ_mean] # x_bar_t |
| 70 | + |
| 71 | + if t == 2 |
| 72 | + first_two = [θ_mean_old'; θ_new'] # θ_mean_old here only contains the first θ |
| 73 | + wum[:vars] = diag(cov(first_two)) |
| 74 | + else#if t <= 1000 |
| 75 | + D = length(θ_new) |
| 76 | + # D = 2.4^2 |
| 77 | + wum[:vars] = (t - 1) / t * wum[:vars] .+ 100 * eps(Float64) + |
| 78 | + (2.4^2 / D) / t * (t * θ_mean_old .* θ_mean_old - (t + 1) * θ_mean_new .* θ_mean_new + θ_new .* θ_new) |
| 79 | + end |
| 80 | + |
| 81 | + if t > 100 |
| 82 | + wum[:stds] = sqrt(wum[:vars]) |
| 83 | + wum[:stds] = wum[:stds] / min(wum[:stds]...) |
| 84 | + end |
46 | 85 | end |
| 86 | +end |
| 87 | + |
| 88 | +update_state(wum::WarmUpManager) = begin |
| 89 | + wum.iter_n += 1 # update iteration number |
47 | 90 |
|
48 | | - if t > 500 |
49 | | - spl.info[:stds] = sqrt(spl.info[:θ_vars]) |
50 | | - spl.info[:stds] = spl.info[:stds] / min(spl.info[:stds]...) |
| 91 | + # Update state |
| 92 | + if wum.state == 1 |
| 93 | + if wum.iter_n > 100 |
| 94 | + wum.state = 2 |
| 95 | + end |
| 96 | + elseif wum.state == 2 |
| 97 | + if wum.iter_n > 900 |
| 98 | + wum.state = 3 |
| 99 | + end |
| 100 | + elseif wum.state == 3 |
| 101 | + if wum.iter_n > 1000 |
| 102 | + wum.state = 4 |
| 103 | + end |
| 104 | + elseif wum.state == 4 |
| 105 | + # no more change |
| 106 | + else |
| 107 | + error("[Turing.WarmUpManager] unknown state $(wum.state)") |
51 | 108 | end |
52 | 109 | end |
53 | 110 |
|
54 | | -init_pre_cond_parameters{T<:Hamiltonian}(vi::VarInfo, spl::Sampler{T}) = begin |
55 | | - spl.info[:θ_mean] = realpart(vi[spl]) |
56 | | - spl.info[:θ_num] = 1 |
57 | | - D = length(vi[spl]) |
58 | | - spl.info[:stds] = ones(D) |
59 | | - spl.info[:θ_vars] = nothing |
| 111 | +adapt(wum::WarmUpManager, stats::Float64, θ_new) = begin |
| 112 | + update_state(wum) |
| 113 | + |
| 114 | + # Use Dual Averaging to adapt ϵ |
| 115 | + if wum.state in [1, 2, 3] |
| 116 | + adapt_step_size(wum, stats) |
| 117 | + end |
| 118 | + |
| 119 | + # Update pre-conditioning matrix |
| 120 | + if wum.state == 2 |
| 121 | + update_pre_cond(wum, θ_new) |
| 122 | + end |
60 | 123 | end |
0 commit comments