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,Ub;n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e3, timeout = Inf) X00 = deepcopy(X0) 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,:] = X0t for i=2:n nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb)); DEBUG && @assert !any(nest[i,:] .> Ub) DEBUG && @assert !any(nest[i,:] .< Lb) end # Get the current best fitness=10^20*ones(n,1); fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness); DEBUG && println("f(X0) = \$(f(X00)), f(bestnest) = \$(fmin)") DEBUG && @assert X00 == X0 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 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 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); DEBUG && @assert !any(nest[j,:] .> Ub) DEBUG && @assert !any(nest[j,:] .< Lb) end return 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,:]; return 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.Ub; s[J] = Ub[J]; return s end # # ## You can replace the following by your own functions # # # A d-dimensional objective function # # function fobj(u) # # ## d-dimensional sphere function sum_j=1^d (u_j-1)^2. # # # with a minimum at (1,1, ...., 1); # # z=sum((u-1).^2); # # end # # dims = 10 # # cuckoo_search(fobj,zeros(dims),-10*ones(dims),10*ones(dims))