Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • mattias
  • sobolfast
3 results

cuckooSearch.jl

Blame
  • cuckooSearch.jl 5.40 KiB
    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