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


keys(Optim.UnconstrainedProblems.examples) = AbstractString["Exponential","Fletcher-Powell","Rosenbrock","Parabola","Hosaki","Large Polynomial","Polynomial","Powell","Himmelblau"]
Out[18]:
ezcontour (generic function with 1 method)

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")


WARNING: Method definition mcmc(Function, Array{T<:Any, 1}, Real, Int64) in module Main at In[33]:8 overwritten at In[34]:8.
WARNING: replacing docs for 'mcmc :: Tuple{Function,Array{T,1},Real,Int64}' in module 'Main'.
mcmc!(f,[0.0;1.0],0.15,hist) = (0.0006545841687509579,[0.97749,0.956703])
WARNING: Method definition mcmc!(Function, Array{T<:Any, 1}, Real, Array{T<:Any, 2}) in module Main at In[33]:12 overwritten at In[34]:12.
WARNING: Method definition mcmc!(Function, Array{T<:Any, 1}, Real, Int64, Array{T<:Any, 2}) in module Main at In[33]:16 overwritten at In[34]:16.
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")


WARNING: Method definition anneal(Function, Array{T<:Any, 1}, Real, Any) in module Main at In[51]:8 overwritten at In[52]:8.
WARNING: replacing docs for 'anneal :: Tuple{Function,Array{T,1},Real,Any}' in module 'Main'.
anneal!(f,[0.0;1.0],0.01,[linspace(Tmax,0,N0);zeros(N1)],hist) = (1.9908015474582628e-7,[1.00028,1.00052])
WARNING: Method definition anneal!(Function, Array{T<:Any, 1}, Real, Any, Array{T<:Any, 2}) in module Main at In[51]:13 overwritten at In[52]:13.
Out[52]:

In [53]:
plot(hist[1,:])


Out[53]: