cuckooSearch.jl 6.18 KB
Newer Older
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
"""
18
function cuckoo_search(f,X0, Lb,Ub;n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e3, timeout = Inf)
19
    X00 = deepcopy(X0)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
20
21
22
23
24
25
26
27
28
29
30
31
32
    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)
33
    nest[1,:] = X0t
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
34
35
    for i=2:n
        nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
36
37
        DEBUG && @assert !any(nest[i,:] .> Ub)
        DEBUG && @assert !any(nest[i,:] .< Lb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
38
39
    end
    # Get the current best
40

Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
41
42
    fitness=10^20*ones(n,1);
    fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness);
43
44
    DEBUG && println("f(X0) = $(f(X00)), f(bestnest) = $(fmin)")
    DEBUG && @assert X00 == X0
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
    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);
81
    println("f(bestnest) = $(fmin)")
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    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);
116
117
        DEBUG && @assert !any(nest[j,:] .> Ub)
        DEBUG && @assert !any(nest[j,:] .< Lb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
118
    end
119
    return nest
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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,:];
135
    return fmin,best,nest,fitness
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

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;
152
153
154
155
    for j = 1:size(nest,1)
        new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub);
    end
    return new_nest
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
156
157
158
159
160
end

# Application of simple constraints
function simplebounds(s,Lb,Ub)
    # Apply the lower bound
161
162
    I = s.<Lb;
    s[I] = Lb[I];
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
163
    # Apply the upper bounds
164
165
166
    J = s.>Ub;
    s[J] = Ub[J];
    return s
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
167
168
169
170
171
172
173
174
175
176
177
178
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))