cmaes.jl

module CMAES

    using Distributions

    ∑(x) = sum(x)
    
    function init_constants(xinit,λ,w,μ)
        D = length(xinit)
    
        μ_w = 1.0 / sum(w.^2)
        c_σ =  (μ_w + 2.0) / (D + μ_w + 5.0)
        
        d_σ = 1.0 + c_σ + 2.0*max(0, √((μ_w-1)/(D+1)) -1) 
        
        c_c = (4 + μ_w/D)/(D + 4 + 2μ_w/D) 
        c_1 =  2/((D+1.3)^2 + μ_w)
        c_μ = min(1-c_1,2*(μ_w-2+1/μ_w)/((D+2)^2 +μ_w)) 
        
        D, μ_w, c_σ, d_σ, c_c, c_1, c_μ
    end
    
    function cmaes(f,xinit, Niter, λ=16, σ=1)
        _weights(μ) = [(log(μ+1)-log(i)) for i=1:μ] / ∑( log(μ+1)-log(j) for j=1:μ )
        
            μ = floor(Int,λ/2)
        w = _weights(μ)
        
        D, μ_w, c_σ, d_σ, c_c, c_1, c_μ = init_constants(xinit,λ,w,μ)
        p_σ, p_c = zeros(D), zeros(D)
        C = diagm(ones(D))    
        
        #
        m = xinit
        x = [zeros(D) for i=1:μ]
        
        Normal(C) = rand(MultivariateNormal(zeros(D),C)) 
        chi_D = √D*(1-1/(4*D) + 1/(21*D^2))
        
        fx = zeros(λ)
        
        for t = 1:Niter
                
            x = [m + σ*Normal(C) for i=1:λ]
            fx = map(f,x)
            
            x = x[sortperm(fx)]
            fx = sort(fx)
            
            m_t = m
            m = ∑( w[i] * x[i] for i=1:μ)
            
            Δm = m-m_t
            
            p_σ = (1-c_σ)*p_σ  + √(c_σ*(2-c_σ)*μ_w) * C^(-1/2) * Δm/σ
            
            h_σ = norm(p_σ) < √(1-(1-c_σ)^(2*(t+1)))*(1.4 + 2/(D+1)) ? 1.0 : 0.0
            
            p_c = (1-c_c)*p_c + h_σ*√(c_σ*(2 - c_σ)*μ_w) * Δm/σ
                
            C = (1-c_1-c_µ + (1-h_σ)*c_1*c_c*(2-c_c)) * C +
                c_1 * p_c * p_c' +
                c_μ * ∑( w[i]/σ^2*( (x[i]-m)  * (x[i]-m)') for i=1:μ)
        
            C = (C+C')/2 #keep symmetric part
                
                σ = σ * exp(c_σ/d_σ*(norm(p_σ)/chi_D -1))
                
            norm(p_σ) < 1e-60 && return x[1],fx[1]
        end
            
        x[1],fx[1]
    end

end


using BlackBoxOptimizationBenchmarking

f = BlackBoxOptimizationBenchmarking.F18
xinit = rand(3)

λ=32
σ=2
Niter = 10_000 

xmin, fmin = CMAES.cmaes(f, xinit, Niter, λ, σ) 

fmin < f.f_opt + 1e-6