cuckooSearch.jl

using Devectorize
"""
`cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)`\n
`n` = Number of nests (or different solutions)
`pa=0.25` Discovery rate of alien eggs/solutions
Change this if you want to get better results
Based on implementation by
@inproceedings{yang2009cuckoo,
  title={Cuckoo search via L{\'e}vy flights},
  author={Yang, Xin-She and Deb, Suash},
  booktitle={Nature \& Biologically Inspired Computing, 2009. NaBIC 2009. World Congress on},
  pages={210--214},
  year={2009},
  organization={IEEE}
}
http://www.mathworks.com/matlabcentral/fileexchange/29809-cuckoo-search--cs--algorithm
"""
function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)
    nd=size(X0,1);
    X0t = X0'
    Lb = Lb'
    Ub = Ub'
    if !all(isfinite(Lb))
        Lb=X0t-0.99999*abs(X0t);
    end
    if !all(isfinite(Ub))
        Ub=X0t+0.99999*abs(X0t);
    end

    # Random initial solutions
    nest = zeros(n,nd)
    nest[1,:] = X0
    for i=2:n
        nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
    end
    # Get the current best
    fitness=10^20*ones(n,1);
    fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness);
    N_iter=0;
    t0 = time()
    ## Starting iterations
    while fmin>Tol && N_iter < max_iter
        # Generate new solutions (but keep the current best)
        new_nest=get_cuckoos(nest,bestnest,Lb,Ub);
        fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness);
        # Update the counter
        N_iter += n;
        if fnew<fmin
            fmin=fnew;
            bestnest=best;
        end
        if time()-t0 > timeout
            display("Cuckoo search: timeout $(timeout)s reached ($(time()-t0)s)")
            break
        end
        # Discovery and randomization
        new_nest=empty_nests(nest,Lb,Ub,pa) ;
        # Evaluate this set of solutions
        fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness);
        # Update the counter again
        N_iter += n;
        # Find the best objective so far
        if fnew<fmin
            fmin=fnew;
            bestnest=best;
        end
        if time()-t0 > timeout
            display("Cuckoo search: timeout $(timeout)s reached ($(time()-t0)s)")
            break
        end
    end ## End of iterations
    ## Post-optimization processing
    ## Display all the nests
    println("Total number of iterations=",N_iter);
    println("f(bestnest) = $(fmin)")
    squeeze(bestnest',2),fmin

end
## --------------- All subfunctions are list below ------------------
## Get cuckoos by ramdom walk
function get_cuckoos(nest,best,Lb,Ub)
    # Levy flights
    n=size(nest,1);
    # Levy exponent and coefficient
    # For details, see equation (2.21), Page 16 (chapter 2) of the book
    # X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
    beta=3/2;
    sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
    for j=1:n
        s=nest[j,:];
        # This is a simple way of implementing Levy flights
        # For standard random walks, use step=1;
        ## Levy flights by Mantegna’s algorithm
        u=randn(size(s))*sigma;
        v=randn(size(s));
        betai = 1/beta
        @devec step=u./abs(v).^betai;
        # In the next equation, the difference factor (s-best) means that
        # when the solution is the best solution, it remains unchanged.
        stepsize=0.01*step.*(s-best);
        # Here the factor 0.01 comes from the fact that L/100 should the typical
        # step size of walks/flights where L is the typical lenghtscale;
        # otherwise, Levy flights may become too aggresive/efficient,
        # which makes new solutions (even) jump out side of the design domain
        # (and thus wasting evaluations).
        # Now the actual random walks or flights
        s=s+stepsize.*randn(size(s));
        # Apply simple bounds/limits
        nest[j,:]=simplebounds(s,Lb,Ub);
    end
    nest
end
## Find the current best nest
function get_best_nest(f,nest,newnest,fitness)
    # Evaluating all new solutions
    for j=1:size(nest,1)
        fnew=f(squeeze(newnest[j,:]',2));
        if fnew<=fitness[j]
            fitness[j]=fnew;
            nest[j,:]=newnest[j,:];

        end
    end
    # Find the current best
    (fmin,K) = findmin(fitness) ;
    best=nest[K,:];
    fmin,best,nest,fitness

end

## Replace some nests by constructing new solutions/nests
function empty_nests(nest,Lb,Ub,pa)
    # A fraction of worse nests are discovered with a probability pa
    n=size(nest,1);
    # Discovered or not -- a status vector
    K=rand(size(nest)).>pa;
    # In the real world, if a cuckoo’s egg is very similar to a host’s eggs, then
    # this cuckoo’s egg is less likely to be discovered, thus the fitness should
    # be related to the difference in solutions. Therefore, it is a good idea
    # to do a random walk in a biased way with some random step sizes.
    ## New solution by biased/selective random walks
    stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]);
    new_nest=nest+stepsize.*K;
    for j = 1:size(nest,1)
        new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub);
    end
    return new_nest
end

# Application of simple constraints
function simplebounds(s,Lb,Ub)
    # Apply the lower bound
    I = s.<Lb;
    s[I] = Lb[I];
    # Apply the upper bounds
    J = s.>Ub;
    s[J] = Ub[J];
    return s
end