In [18]:
using Optim
@show keys(Optim.UnconstrainedProblems.examples)
using Plots
plotlyjs()
ezcontour(x, y, f) = begin
X = repmat(x', length(y), 1)
Y = repmat(y, 1, length(x))
# Evaluate each f(x, y)
Z = map((x,y) -> (f([x,y])), X, Y)
plot(x, y, Z, st=:contour)
end
Out[18]:
In [34]:
""" Run a Markov Chain Monte Carlo scheme to sample from
a distribution that minimizes a function f
mcmc!(f, x0, tau, hist) -> save the history in hist, which determins niter
mcmc(f, x0, tau, niter) -> return the best point
"""
function mcmc(f::Function, x0::Vector, tau::Real, niter::Int)
mcmc!(f,x0,tau,niter,zeros(0,0))
end
function mcmc!(f::Function, x0::Vector, tau::Real, hist::Matrix)
mcmc!(f, x0, tau, size(hist,2)-1, hist)
end
function mcmc!(f::Function, x0::Vector, tau::Real, niter::Int, hist::Matrix)
if !isempty(hist)
@assert niter == size(hist,2)-1
@assert length(x0) == size(hist,1)-1
end
n = length(x0)
bestx = x0
bestf = f(x0)
if !isempty(hist)
# save the point
hist[:,1] = [bestf;x0]
end
curf = bestf
x = bestx
for iter=1:niter
xn = x + tau*randn(n)
nextf = f(xn)
a = exp(-nextf)/exp(-curf)
if a > 1 || rand() <= a
x = xn
curf = nextf
end
if curf < bestf
bestx = x
bestf = curf
end
if !isempty(hist)
# save the point
hist[:,iter+1] = [curf; x]
end
end
return bestf, bestx
end
f = Optim.UnconstrainedProblems.examples["Rosenbrock"].f
hist = zeros(3,10000)
@show mcmc!(f, [0.0;1.0], 0.15, hist)
ezcontour(-1:0.05:2, -1:0.05:4, x -> (log10(f(x))))
plot!(hist[2,:],hist[3,:],label="",color="darkblue")
Out[34]:
In [52]:
""" Run a Markov Chain Monte Carlo scheme to sample from
a distribution that minimizes a function f
mcmc!(f, x0, tau, hist) -> save the history in hist, which determins niter
mcmc(f, x0, tau, niter) -> return the best point
"""
function anneal(f::Function, x0::Vector, tau::Real, tvals)
mcmc!(f,x0,tau,tvals,zeros(0,0))
end
function anneal!(f::Function, x0::Vector, tau::Real, tvals, hist::Matrix)
niter = length(tvals)
if !isempty(hist)
@assert niter == size(hist,2)-1
@assert length(x0) == size(hist,1)-1
end
n = length(x0)
bestx = x0
bestf = f(x0)
if !isempty(hist)
# save the point
hist[:,1] = [bestf;x0]
end
curf = bestf
x = bestx
for iter=1:niter
xn = x + tau*randn(n)
nextf = f(xn)
a = exp((-nextf + curf)/tvals[iter])
#@show a
if a > 1 || rand() <= a
x = xn
curf = nextf
end
if curf < bestf
bestx = x
bestf = curf
end
if !isempty(hist)
# save the point
hist[:,iter+1] = [curf; x]
end
end
return bestf, bestx
end
f = Optim.UnconstrainedProblems.examples["Rosenbrock"].f
Tmax = 1
N0 = 5000
N1 = 5000
hist = zeros(3,N0+N1+1)
@show anneal!(f, [0.0;1.0], 0.01, [linspace(Tmax,0,N0); zeros(N1)], hist)
ezcontour(-1:0.05:2, -1:0.05:4, x -> (log10(f(x))))
plot!(hist[2,:],hist[3,:],label="",color="darkblue")
Out[52]:
In [53]:
plot(hist[1,:])
Out[53]: