diff --git a/src/cuckooSearch.jl b/src/cuckooSearch.jl
new file mode 100644
index 0000000000000000000000000000000000000000..0fed3056a9fca939415471597755ac0a7b3b60f2
--- /dev/null
+++ b/src/cuckooSearch.jl
@@ -0,0 +1,167 @@
+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);
+ 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;
+end
+
+# Application of simple constraints
+function simplebounds(s,Lb,Ub)
+ # Apply the lower bound
+ ns_tmp=s;
+ I=ns_tmp.<Lb;
+ ns_tmp[I]=Lb[I];
+ # Apply the upper bounds
+ J=ns_tmp.>Ub;
+ ns_tmp[J]=Ub[J];
+ # Update this new move
+ s=ns_tmp;
+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))