بهینه سازی (ریاضیاتی) چیست؟ — راهنمای جامع

1def euc_2d(c1, c2)
2  Math.sqrt((c1[0] - c2[0])**2.0 + (c1[1] - c2[1])**2.0).round
3end
4
5def cost(permutation, cities)
6  distance =0
7  permutation.each_with_index do |c1, i|
8    c2 = (i==permutation.size-1) ? permutation[0] : permutation[i+1]
9    distance += euc_2d(cities[c1], cities[c2])
10  end
11  return distance
12end
13
14def random_permutation(cities)
15  perm = Array.new(cities.size){|i| i}
16  perm.each_index do |i|
17    r = rand(perm.size-i) + i
18    perm[r], perm[i] = perm[i], perm[r]
19  end
20  return perm
21end
22
23def calculate_neighbor_rank(city_number, cities, ignore=[])
24  neighbors = []
25  cities.each_with_index do |city, i|
26    next if i==city_number or ignore.include?(i)
27    neighbor = {:number=>i}
28    neighbor[:distance] = euc_2d(cities[city_number], city)
29    neighbors << neighbor
30  end
31  return neighbors.sort!{|x,y| x[:distance] <=> y[:distance]}
32end
33
34def get_edges_for_city(city_number, permutation)
35  c1, c2 = nil, nil
36  permutation.each_with_index do |c, i|
37    if c == city_number
38      c1 = (i==0) ? permutation.last : permutation[i-1]
39      c2 = (i==permutation.size-1) ? permutation.first : permutation[i+1]
40      break
41    end
42  end
43  return [c1, c2]
44end
45
46def calculate_city_fitness(permutation, city_number, cities)
47  c1, c2 = get_edges_for_city(city_number, permutation)
48  neighbors = calculate_neighbor_rank(city_number, cities)
49  n1, n2 = -1, -1
50  neighbors.each_with_index do |neighbor,i|
51    n1 = i+1 if neighbor[:number] == c1
52    n2 = i+1 if neighbor[:number] == c2
53    break if n1!=-1 and n2!=-1
54  end
55  return 3.0 / (n1.to_f + n2.to_f)
56end
57
58def calculate_city_fitnesses(cities, permutation)
59  city_fitnesses = []
60  cities.each_with_index do |city, i|
61    city_fitness = {:number=>i}
62    city_fitness[:fitness] = calculate_city_fitness(permutation, i, cities)
63    city_fitnesses << city_fitness
64  end
65  return city_fitnesses.sort!{|x,y| y[:fitness] <=> x[:fitness]}
66end
67
68def calculate_component_probabilities(ordered_components, tau)
69  sum = 0.0
70  ordered_components.each_with_index do |component, i|
71    component[:prob] = (i+1.0)**(-tau)
72    sum += component[:prob]
73  end
74  return sum
75end
76
77def make_selection(components, sum_probability)
78  selection = rand()
79  components.each_with_index do |component, i|
80    selection -= (component[:prob] / sum_probability)
81    return component[:number] if selection <= 0.0
82  end
83  return components.last[:number]
84end
85
86def probabilistic_selection(ordered_components, tau, exclude=[])
87  sum = calculate_component_probabilities(ordered_components, tau)
88  selected_city = nil
89  begin
90    selected_city = make_selection(ordered_components, sum)
91  end while exclude.include?(selected_city)
92  return selected_city
93end
94
95def vary_permutation(permutation, selected, new, long_edge)
96  perm = Array.new(permutation)
97  c1, c2 = perm.rindex(selected), perm.rindex(new)
98  p1,p2 = (c1<c2) ? [c1,c2] : [c2,c1]
99  right = (c1==perm.size-1) ? 0 : c1+1
100  if perm[right] == long_edge
101    perm[p1+1..p2] = perm[p1+1..p2].reverse
102  else
103    perm[p1...p2] = perm[p1...p2].reverse
104  end
105  return perm
106end
107
108def get_long_edge(edges, neighbor_distances)
109  n1 = neighbor_distances.find {|x| x[:number]==edges[0]}
110  n2 = neighbor_distances.find {|x| x[:number]==edges[1]}
111  return (n1[:distance] > n2[:distance]) ? n1[:number] : n2[:number]
112end
113
114def create_new_perm(cities, tau, perm)
115  city_fitnesses = calculate_city_fitnesses(cities, perm)
116  selected_city = probabilistic_selection(city_fitnesses.reverse, tau)
117  edges = get_edges_for_city(selected_city, perm)
118  neighbors = calculate_neighbor_rank(selected_city, cities)
119  new_neighbor = probabilistic_selection(neighbors, tau, edges)
120  long_edge = get_long_edge(edges, neighbors)
121  return vary_permutation(perm, selected_city, new_neighbor, long_edge)
122end
123
124def search(cities, max_iterations, tau)
125  current = {:vector=>random_permutation(cities)}
126  current[:cost] = cost(current[:vector], cities)
127  best = current
128  max_iterations.times do |iter|
129    candidate = {}
130    candidate[:vector] = create_new_perm(cities, tau, current[:vector])
131    candidate[:cost] = cost(candidate[:vector], cities)
132    current = candidate
133    best = candidate if candidate[:cost] < best[:cost]
134    puts " > iter #{(iter+1)}, curr=#{current[:cost]}, best=#{best[:cost]}"
135  end
136  return best
137end
138
139if __FILE__ == $0
140  # problem configuration
141  berlin52 = [[565,575],[25,185],[345,750],[945,685],[845,655],
142   [880,660],[25,230],[525,1000],[580,1175],[650,1130],[1605,620],
143   [1220,580],[1465,200],[1530,5],[845,680],[725,370],[145,665],
144   [415,635],[510,875],[560,365],[300,465],[520,585],[480,415],
145   [835,625],[975,580],[1215,245],[1320,315],[1250,400],[660,180],
146   [410,250],[420,555],[575,665],[1150,1160],[700,580],[685,595],
147   [685,610],[770,610],[795,645],[720,635],[760,650],[475,960],
148   [95,260],[875,920],[700,500],[555,815],[830,485],[1170,65],
149   [830,610],[605,625],[595,360],[1340,725],[1740,245]]
150  # algorithm configuration
151  max_iterations = 250
152  tau = 1.8
153  # execute the algorithm
154  best = search(berlin52, max_iterations, tau)
155  puts "Done. Best Solution: c=#{best[:cost]}, v=#{best[:vector].inspect}"
156end

1def onemax(vector)
2  return vector.inject(0){|sum, value| sum + value}
3end
4
5def random_bitstring(size)
6  return Array.new(size){ ((rand()<0.5) ? 1 : 0) }
7end
8
9def path_exists?(i, j, graph)
10  visited, stack = [], [i]
11  while !stack.empty?
12    return true if stack.include?(j)
13    k = stack.shift
14    next if visited.include?(k)
15    visited << k
16    graph[k][:out].each {|m| stack.unshift(m) if !visited.include?(m)}
17  end
18  return false
19end
20
21def can_add_edge?(i, j, graph)
22  return !graph[i][:out].include?(j) && !path_exists?(j, i, graph)
23end
24
25def get_viable_parents(node, graph)
26  viable = []
27  graph.size.times do |i|
28    if node!=i and can_add_edge?(node, i, graph)
29      viable << i
30    end
31  end
32  return viable
33end
34
35def compute_count_for_edges(pop, indexes)
36  counts = Array.new(2**(indexes.size)){0}
37  pop.each do |p|
38    index = 0
39    indexes.reverse.each_with_index do |v,i|
40      index += ((p[:bitstring][v] == 1) ? 1 : 0) * (2**i)
41    end
42    counts[index] += 1
43  end
44 return counts
45end
46
47def fact(v)
48  return v <= 1 ? 1 : v*fact(v-1)
49end
50
51def k2equation(node, candidates, pop)
52  counts = compute_count_for_edges(pop, [node]+candidates)
53  total = nil
54  (counts.size/2).times do |i|
55    a1, a2 = counts[i*2], counts[(i*2)+1]
56    rs = (1.0/fact((a1+a2)+1).to_f) * fact(a1).to_f * fact(a2).to_f
57    total = (total.nil? ? rs : total*rs)
58  end
59  return total
60end
61
62def compute_gains(node, graph, pop, max=2)
63  viable = get_viable_parents(node[:num], graph)
64  gains = Array.new(graph.size) {-1.0}
65  gains.each_index do |i|
66    if graph[i][:in].size < max and viable.include?(i)
67      gains[i] = k2equation(node[:num], node[:in]+[i], pop)
68    end
69  end
70  return gains
71end
72
73def construct_network(pop, prob_size, max_edges=3*pop.size)
74  graph = Array.new(prob_size) {|i| {:out=>[], :in=>[], :num=>i} }
75  gains = Array.new(prob_size)
76  max_edges.times do
77    max, from, to = -1, nil, nil
78    graph.each_with_index do |node, i|
79      gains[i] = compute_gains(node, graph, pop)
80      gains[i].each_with_index {|v,j| from,to,max = i,j,v if v>max}
81    end
82    break if max <= 0.0
83    graph[from][:out] << to
84    graph[to][:in] << from
85  end
86  return graph
87end
88
89def topological_ordering(graph)
90  graph.each {|n| n[:count] = n[:in].size}
91  ordered,stack = [], graph.select {|n| n[:count]==0}
92  while ordered.size < graph.size
93    current = stack.shift
94    current[:out].each do |edge|
95      node = graph.find {|n| n[:num]==edge}
96      node[:count] -= 1
97      stack << node if node[:count] <= 0
98    end
99    ordered << current
100  end
101  return ordered
102end
103
104def marginal_probability(i, pop)
105  return pop.inject(0.0){|s,x| s + x[:bitstring][i]} / pop.size.to_f
106end
107
108def calculate_probability(node, bitstring, graph, pop)
109  return marginal_probability(node[:num], pop) if node[:in].empty?
110  counts = compute_count_for_edges(pop, [node[:num]]+node[:in])
111  index = 0
112  node[:in].reverse.each_with_index do |v,i|
113    index += ((bitstring[v] == 1) ? 1 : 0) * (2**i)
114  end
115  i1 = index + (1*2**(node[:in].size))
116  i2 = index + (0*2**(node[:in].size))
117  a1, a2 = counts[i1].to_f, counts[i2].to_f
118  return a1/(a1+a2)
119end
120
121def probabilistic_logic_sample(graph, pop)
122  bitstring = Array.new(graph.size)
123  graph.each do |node|
124    prob = calculate_probability(node, bitstring, graph, pop)
125    bitstring[node[:num]] = ((rand() < prob) ? 1 : 0)
126  end
127  return {:bitstring=>bitstring}
128end
129
130def sample_from_network(pop, graph, num_samples)
131  ordered = topological_ordering(graph)
132  samples = Array.new(num_samples) do
133    probabilistic_logic_sample(ordered, pop)
134  end
135  return samples
136end
137
138def search(num_bits, max_iter, pop_size, select_size, num_children)
139  pop = Array.new(pop_size) { {:bitstring=>random_bitstring(num_bits)} }
140  pop.each{|c| c[:cost] = onemax(c[:bitstring])}
141  best = pop.sort!{|x,y| y[:cost] <=> x[:cost]}.first
142  max_iter.times do |it|
143    selected = pop.first(select_size)
144    network = construct_network(selected, num_bits)
145    arcs = network.inject(0){|s,x| s+x[:out].size}
146    children = sample_from_network(selected, network, num_children)
147    children.each{|c| c[:cost] = onemax(c[:bitstring])}
148    children.each {|c| puts " >>sample, f=#{c[:cost]} #{c[:bitstring]}"}
149    pop = pop[0...(pop_size-select_size)] + children
150    pop.sort! {|x,y| y[:cost] <=> x[:cost]}
151    best = pop.first if pop.first[:cost] >= best[:cost]
152    puts " >it=#{it}, arcs=#{arcs}, f=#{best[:cost]}, [#{best[:bitstring]}]"
153    converged = pop.select {|x| x[:bitstring]!=pop.first[:bitstring]}.empty?
154    break if converged or best[:cost]==num_bits
155  end
156  return best
157end
158
159if __FILE__ == $0
160  # problem configuration
161  num_bits = 20
162  # algorithm configuration
163  max_iter = 100
164  pop_size = 50
165  select_size = 15
166  num_children = 25
167  # execute the algorithm
168  best = search(num_bits, max_iter, pop_size, select_size, num_children)
169  puts "done! Solution: f=#{best[:cost]}/#{num_bits}, s=#{best[:bitstring]}"
170end

1clc, clear, close all
2%% Set problem parameters
3% select your cost function:
4costFunction = 'rastrigin';           % F6 +/-5.12
5npar  = 100;          % number of optimization variables
6varLo = -5;         % Lower  band of parameter
7varHi =  5;         % Higher band of parameter
8%% Set COA parameters
9numCuckooS = 5;             % number of initial population
10minNumberOfEggs = 2;        % minimum number of eggs for each cuckoo
11maxNumberOfEggs = 4;        % maximum number of eggs for each cuckoo
12maxIter = 100;              % maximum iterations of the Cuckoo Algorithm
13knnClusterNum = 1;          % number of clusters that we want to make
14motionCoeff = 9;            % Lambda variable in COA paper, default=2
15accuracy = -inf;           % How much accuracy in answer is needed
16maxNumOfCuckoos = 10;      % maximum number of cuckoos that can live at the same time
17radiusCoeff = 5;           % Control parameter of egg laying
18cuckooPopVariance = 1e-13;   % population variance that cuts the optimization
19%% initialize population:
20cuckooPop = cell(numCuckooS,1);
21% initialize egg laying center for each cuckoo
22for cuckooNumber = 1:numCuckooS    
23    cuckooPop{cuckooNumber}.center = ( varHi-varLo )*rand(1,npar) + varLo;
24end
25%% initialize COA cost minimization plot
26figure(1)
27hold on
28xlabel 'Cuckoo iteration'
29ylabel 'Cost Value'
30%% Start Cuckoo Optimization Algorithm
31iteration = 0;
32maxProfit = -1e20;        % Let initial value be negative number
33goalPoint = (varHi - varLo)*rand(1,npar) + varLo; % a random goalpoint to start COA
34globalBestCuckoo = goalPoint;
35globalMaxProfit = maxProfit;
36profitVector = [];
37while ( (iteration <= maxIter) && (-maxProfit > accuracy) )
38    
39    iteration = iteration + 1
40    
41    % initialize number of eggs for each cuckoo
42    for cuckooNumber = 1:numCuckooS        
43        cuckooPop{cuckooNumber}.numberOfEggs = floor( (maxNumberOfEggs - minNumberOfEggs) * rand + minNumberOfEggs );
44    end
45    % get total number of available eggs between all cuckooS
46    summ = 0;
47    for cuckooNumber = 1:numCuckooS
48        summ = summ + cuckooPop{cuckooNumber}.numberOfEggs;
49    end
50    % calculate egg laying radius for each Cuckoo, considering problem
51    % limitations and ratio of each cuckoo's eggs
52    for cuckooNumber = 1:numCuckooS
53        cuckooPop{cuckooNumber}.eggLayingRadius = cuckooPop{cuckooNumber}.numberOfEggs/summ * ( radiusCoeff * (varHi-varLo) );
54    end
55    % To lay eggs, we produced some radius values less than egg laying
56    % radius of each cuckoo
57    for cuckooNumber = 1:numCuckooS
58        cuckooPop{cuckooNumber}.eggLayingRadiuses = cuckooPop{cuckooNumber}.eggLayingRadius * rand(cuckooPop{cuckooNumber}.numberOfEggs,1);
59    end
60    
61    for cuckooNumber = 1:numCuckooS
62        params = cuckooPop{cuckooNumber}.center;        % get center values
63        tmpRadiuses = cuckooPop{cuckooNumber}.eggLayingRadiuses;
64        numRadiuses = numel(tmpRadiuses);
65        % divide a (hyper)circle to 'numRadiuses' segments
66        % This is to search all over the current habitat
67        angles = linspace(0,2*pi,numRadiuses);    % in Radians
68        newParams = [];
69        for cnt = 1:numRadiuses
70            addingValue = zeros(1,npar);
71            for iii = 1:npar
72                randNum = floor(2*rand)+1;
73                addingValue(iii) = (-1)^randNum * tmpRadiuses(cnt)*cos(angles(cnt)) + tmpRadiuses(cnt)*sin(angles(cnt));
74            end
75            newParams = [newParams; params +  addingValue ];
76        end
77        
78        
79        % check for variable limits
80        newParams(find(newParams>varHi)) = varHi;
81        newParams(find(newParams<varLo)) = varLo;
82        cuckooPop{cuckooNumber}.newPosition4Egg = newParams;
83    end
84    
85    % Now egg laying is done
86    
87    % Now that egg positions are found, they are laid, and so its time to
88    % remove the eggs on the same positions (because each egg only can go to one nest)
89    for cuckooNumber = 1:numCuckooS
90        tmpPopulation = cuckooPop{cuckooNumber}.newPosition4Egg;
91        tmpPopulation = floor(tmpPopulation * 1e20)/1e20;
92        ii = 2;
93        cntt = 1;
94        while ii <= size(tmpPopulation,1) || cntt <= size(tmpPopulation,1)
95            if all((tmpPopulation(cntt,:) == tmpPopulation(ii,:)))
96                tmpPopulation(ii,:) = [];
97            end
98            ii = ii + 1;
99            if ii > size(tmpPopulation,1) && cntt <= size(tmpPopulation,1)
100                cntt = cntt + 1;
101                ii = cntt + 1;
102                if ii > size(tmpPopulation,1)
103                    break
104                end
105            end
106        end
107        cuckooPop{cuckooNumber}.newPosition4Egg = tmpPopulation;
108    end    
109    
110     
111    % Now we evalute egg positions
112    for cuckooNumber = 1:numCuckooS
113        cuckooPop{cuckooNumber}.profitValues = -feval(costFunction,[cuckooPop{cuckooNumber}.center ; cuckooPop{cuckooNumber}.newPosition4Egg]);        
114    end
115    
116    % Now we check to see if cuckoo population is more than maxNumOfCuckoos
117    % this case we should keep 1st maxNumOfCuckoos cuckoos and kill the others
118    allPositions = [];
119    whichCuckooPopTheEggBelongs = [];
120    tmpProfits = [];
121    if numCuckooS > maxNumOfCuckoos
122        for cuckooNumber = 1:numCuckooS
123            tmpProfits = [tmpProfits; cuckooPop{cuckooNumber}.profitValues];
124            allPositions = [allPositions; [cuckooPop{cuckooNumber}.center; cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar)]];
125            whichCuckooPopTheEggBelongs = [whichCuckooPopTheEggBelongs; cuckooNumber*ones(size(cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar),1),1) ];
126        end
127        % now we sort cuckoo profits
128        [sortedProfits, sortingIndex] = sort(tmpProfits,'descend');
129        % Get best cuckoo to be copied to next generation
130        bestCuckooCenter = allPositions(sortingIndex(1),1:npar);
131        
132        sortedProfits = sortedProfits(1:maxNumOfCuckoos);
133        allPositions = allPositions(sortingIndex(1:maxNumOfCuckoos),:);
134        clear cuckooPop
135        for ii = 1:maxNumOfCuckoos
136            cuckooPop{ii}.newPosition4Egg = allPositions(ii,:);
137            cuckooPop{ii}.center = allPositions(ii,:);
138            cuckooPop{ii}.profitValues = sortedProfits(ii);
139        end
140        numCuckooS = maxNumOfCuckoos;
141    else
142        for cuckooNumber = 1:numCuckooS
143            tmpProfits = [tmpProfits; cuckooPop{cuckooNumber}.profitValues];
144            allPositions = [allPositions; [cuckooPop{cuckooNumber}.center; cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar)] ];
145            whichCuckooPopTheEggBelongs = [whichCuckooPopTheEggBelongs; cuckooNumber*ones(size(cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar),1),1) ];
146        end
147        [sortedProfits, sortingIndex] = sort(tmpProfits,'descend');
148        % Get best cuckoo to be copied to next generation
149        bestCuckooCenter = allPositions(sortingIndex(1),1:npar);
150    end
151    
152    maxProfit  = sortedProfits(1);
153    currentBestCuckoo = bestCuckooCenter;
154    currentMaxProfit = -feval(costFunction,currentBestCuckoo);
155    if currentMaxProfit > globalMaxProfit
156        globalBestCuckoo = currentBestCuckoo;
157        globalMaxProfit = currentMaxProfit;
158    end
159    
160    % Update cost minimization plot
161    plot(iteration, -globalMaxProfit,'ks','linewidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10)
162    title([ 'Curent Cost = ' num2str(-globalMaxProfit) ' , at Iteration = ' num2str(iteration) ])
163    pause(0.01)
164    
165    profitVector = [profitVector globalMaxProfit];
166    
167    % ======== now we have some eggs that are safe and will grow up ==========
168    %========= mating: =============
169    % first we put all egg positions under each other
170    allPositions = [];
171    whichCuckooPopTheEggBelongs = [];
172    for cuckooNumber = 1:numCuckooS
173        allPositions = [allPositions; cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar)];
174        whichCuckooPopTheEggBelongs = [whichCuckooPopTheEggBelongs; cuckooNumber*ones(size(cuckooPop{cuckooNumber}.newPosition4Egg(:,1:npar),1),1) ];
175    end
176    if sum(var(allPositions)) < cuckooPopVariance
177        break
178    else
179        [clusterNumbers, clusterCenters] = kmeans(allPositions,knnClusterNum);
180    end
181    % make newly made clusters
182    cluster = cell(knnClusterNum,1);
183    for ii = 1:knnClusterNum
184        cluster{ii}.positions = [];
185        cluster{ii}.profits = [];
186    end
187    
188    pointer = zeros(numel(cuckooPop),1);
189    for tmpCounter = 1:numel(cuckooPop)
190        nEggs = size(cuckooPop{tmpCounter}.newPosition4Egg,1);
191        pointer(tmpCounter) = nEggs;
192    end
193    for cnt = 1:length(clusterNumbers)
194        cluster{clusterNumbers(cnt)}.positions = [cluster{clusterNumbers(cnt)}.positions; cuckooPop{whichCuckooPopTheEggBelongs(cnt)}.newPosition4Egg(end-pointer(whichCuckooPopTheEggBelongs(cnt))+1,1:npar)];
195        cluster{clusterNumbers(cnt)}.profits   = [cluster{clusterNumbers(cnt)}.profits; cuckooPop{whichCuckooPopTheEggBelongs(cnt)}.profitValues(end-pointer(whichCuckooPopTheEggBelongs(cnt))+1)];
196        pointer(whichCuckooPopTheEggBelongs(cnt)) = pointer(whichCuckooPopTheEggBelongs(cnt)) - 1;
197    end
198    
199    % Determine the best cluster
200    f_mean = zeros(knnClusterNum,1);
201    for cnt = 1:knnClusterNum
202        f_mean(cnt) = mean(cluster{cnt}.profits);
203    end
204    [sorted_f_mean, sortingIndex_f_mean] = sort(f_mean,'descend');
205    maxFmean = sorted_f_mean(1);   indexOfBestCluster = sortingIndex_f_mean(1);
206    
207    % now that we know group with number 'indexOfBestCluster' is the best we 
208    % should select their best point az Goal Point of other groups
209    [maxProfitInBestCluster, indexOfBestEggPosition] = max(cluster{indexOfBestCluster}.profits);
210    goalPoint  = cluster{indexOfBestCluster}.positions(indexOfBestEggPosition,1:npar);
211    
212    % now all other mature Cuckoos must go toward this goal point for laying
213    % their eggs
214    numNewCuckooS = 0;
215    for cntClstr = 1:size(cluster,1)
216        tmpCluster = cluster{cntClstr};
217        tmpPositions = tmpCluster.positions;
218        for cntPosition = 1:size(tmpPositions,1)
219            tmpPositions(cntPosition,:) = tmpPositions(cntPosition,:) + ...
220                                          motionCoeff * rand(1,npar) .*  (goalPoint  - tmpPositions(cntPosition,:));
221        end
222        % check if variables are in range
223        tmpPositions(find( tmpPositions>varHi )) = varHi;
224        tmpPositions(find( tmpPositions<varLo )) = varLo;
225        % update cluster positions
226        cluster{cntClstr}.positions = tmpPositions;
227        cluster{cntClstr}.center = mean(tmpPositions);
228        % update number of cuckoos: numCuckooS
229        numNewCuckooS = numNewCuckooS + size(cluster{cntClstr}.positions,1);
230    end
231    numCuckooS = numNewCuckooS;
232    % update cuckooPop
233    clear cuckooPop
234    cuckooPop = cell(numCuckooS,1);
235    cntNumCuckooS = 1;
236    for cnt = 1:size(cluster,1)
237        tmpCluster = cluster{cnt};
238        tmpPositions = tmpCluster.positions;
239        for cntPosition = 1:size(tmpPositions,1)
240            cuckooPop{cntNumCuckooS}.center = cluster{cnt}.positions(cntPosition,1:npar);
241            cntNumCuckooS = cntNumCuckooS + 1;
242        end
243    end
244    % Copy the Best cuckoo and its randomized form of this population to go 
245    % to the next generation
246    currentBestCuckoo = bestCuckooCenter;
247    currentMaxProfit = -feval(costFunction,currentBestCuckoo);
248    if currentMaxProfit > globalMaxProfit
249        globalBestCuckoo = currentBestCuckoo;
250        globalMaxProfit = currentMaxProfit;
251    end
252    cuckooPop{end}.center = globalBestCuckoo; % This is because the best cuckoo will live longer and won't die right after egg laying
253    cuckooPop{end}.profitValues = -feval(costFunction,cuckooPop{end}.center);
254    tmp = rand(1,npar).*globalBestCuckoo;
255    tmp(find( tmp>varHi )) = varHi;
256    tmp(find( tmp<varLo )) = varLo;
257    cuckooPop{end-1}.center = tmp;
258    cuckooPop{end-1}.profitValues = -feval(costFunction,cuckooPop{end-1}.center);
259    
260end     % end of while
261%% Now Algorithm has stopped
262disp('Best Params = ')
263disp(globalBestCuckoo')
264fprintf('\n')
265disp('Cost = ')
266disp(-globalMaxProfit)
267% profit history is in:  profitVector
268costVector = - profitVector;
269figure
270plot(costVector,'-ks','linewidth',3,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',15)
271xlabel 'Cuckoo Iteration'
272ylabel 'Cost Value'
273title(['Current Cost = ' num2str(costVector(end)) ', at iteration = ' num2str(iteration) ])

1function y = rastrigin(x)
2y = 10.0 * size(x,2) + sum(x .^2 - 10.0 * cos(2 * pi .* x),2);