Cuckoo search project

.gitignore

@@ -26,3 +26,6 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+
+# Directory generated by VSCode
+.vscode/

Project.toml

@@ -6,4 +6,5 @@ version = "0.1.0"
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/make.jl

@@ -2,12 +2,14 @@ using Documenter
 
 using STMOZOO
 using STMOZOO.Example
+using STMOZOO.Cuckoo
 
 makedocs(sitename="STMO ZOO",
     format = Documenter.HTML(),
-    modules=[Example], # add your module
+    modules=[Example, Cuckoo],  
     pages=Any[
-        "Example"=> "man/example.md",  # add the page to your documentation
+        "Example"=> "man/example.md",  
+        "Cuckoo"=> "man/cuckoo.md"
     ])
 
 #=

docs/src/man/cuckoo.md

@@ -0,0 +1,18 @@
+## Cuckoo search
+
+This is the documentation for the cuckoo search method, implemented by Claudia Mengoni. The module allows to solve a minimization problem through the use of the biological-inspired metaheuristic of cuckoo search.
+
+The method is based on the parasitism of some cuckoo species that exploit the resources of other bird species by laying eggs into their nests.
+To abstract this concept to a computational method the phenomenon is simplified by three main rules:
+* Each cuckoo lays one egg at a time into a randomly chosen nest
+* The nests containing the best eggs are carried over (elitist selection) to the next generation
+* The number of available host nests is fixed and at each generation the alien cuckoo egg can be discovered with a certain probability. At this point the nest is abandoned and a new nest is generated.
+
+It's important to note that this implementation allows a single egg for each nest, hence the terms nest and egg are used interchangeably.
+
+The module contains a function to solve a minimization problem and a helper function that generates an initial population to use as input to the main method.
+
+```@docs 
+init_nests
+cuckoo! 
+```

instructions.md

@@ -71,18 +71,20 @@ Finally, you have to add a [Pluto](https://github.com/fonsp/Pluto.jl) notebook t
 
 Each of you will have to perform a code review of two other projects on **January 7**. You have the full day to do this, though it should not take too long. The aim is to **help** the other groups to make each other's project even better.
 
-- [ ] make a fork of the repo of the person you are reviewing;
+- [ ] clone the projects you have to review, so you can run them locally;
 - [ ] check the source code, is the documentation clear? Anything obvious that can be improved.
 - [ ] run the tests. Do they work? Anything that could be tested but is not done so?
 - [ ] Is the documentation page clear? Do you find any typos? Could an example be added?
 - [ ] Take a look at the notebook. Any suggestions there to improve this?
 
-Big things can be addressed by opening an issue. Small fixes and suggestions to the other person's code can be done immediately and via a pull request.
+You can leave comments in the pull request page. If you have things you want to fix yourself you can also perform a pull request like I did in your code review. **Communicate to the person when you are finised.**
 
 on **January 8**, you have the full day to:
 - [ ] merge the entire request and fix any issues you find meaningful.
 - [ ] fill in a small questionnaire on Ufora about your project and the projects you have reviewed.
 
+Please address all the comments you get as a common courtesy. This does not mean you have to implement all suggestions, but you do have to clarify on why you won't/can't do them.
+
 Afterward, I will approve all pull requests, finalizing the package.
 
 > Two days are blocked for the code review because some of you are in different time zones. This does **not** mean you will need these days in full. You will need approximately reserve one hour for each day. For the remainder, you can do whatever you want.

notebook/cuckoo.jl

@@ -0,0 +1,291 @@
+### A Pluto.jl notebook ###
+# v0.12.11
+
+using Markdown
+using InteractiveUtils
+
+# ╔═╡ 694cf0e0-405a-11eb-08f9-bd9a2e0c6232
+using STMOZOO.Cuckoo, Plots
+
+# ╔═╡ d03cf1ca-4059-11eb-3339-99c511eba3c8
+md"""# Cuckoo search method"""
+
+# ╔═╡ db128a7a-4085-11eb-08b9-5b5199efecc3
+md"""This notebook is meant to show the functioning of the cuckoo search method to solve a minimization problem.
+The method is based on the parasitism of some cuckoo species that exploit the resources of other bird species by laying eggs into their nests.
+
+To abstract this concept to a computational method the phenomenon is simplified into three main rules:
+* Each cuckoo lays one egg at a time into a randomly chosen nest
+* The nests containing the best eggs are carried over (elitist selection) to the next generation
+* The number of available host nests is fixed and at each generation the alien cuckoo egg can be discovered with a certain probability. At this point the nest is abandoned and a new nest is generated.
+
+The peculiarity of this method is that in the first step, the position of each new nest is chosen randomly from a Lévy distribution, mimicking the Lévy flight which is a fliying pattern observed in some bird species.
+
+It's important to note that this implementation allows a single egg for each nest, hence the terms nest and egg are used interchangeably.
+
+More details about the implementation are found in the documentation of the module `STMOZOO.Cuckoo`.
+"""
+
+# ╔═╡ 3c8f220e-4086-11eb-175b-f92b3e3e8253
+md"""## A simple example"""
+
+# ╔═╡ 56c0a760-4086-11eb-1c5f-45828eb40989
+md"""To show the functioning of the method we use the Ackley function in its 2-dimensional form which has its global minimum at (0,0)."""
+
+# ╔═╡ be3ef176-4086-11eb-365c-e1747053c12d
+begin
+	function ackley(x; a=20, b=0.2, c=2π) 
+		d = length(x) 
+		return -a * exp(-b*sqrt(sum(x.^2)/d)) -exp(sum(cos.(c .* x))/d) 
+	end 
+
+	ackley(x...; kwargs...) = ackley(x; kwargs...)
+end
+
+# ╔═╡ e9040d5a-4087-11eb-03a6-51e941460bc5
+md"""This is the landscape of the function: """
+
+# ╔═╡ e0165c02-405a-11eb-296f-05e437e873e3
+pobj = heatmap(-10:0.01:10, -10:0.01:10, ackley)
+
+# ╔═╡ 2607802e-4088-11eb-1a50-bbabd9cbeaee
+md"""We specify a lower and an upper limit for each dimension of the problem we want to solve."""
+
+# ╔═╡ 6cfa745e-408c-11eb-345c-29b71c84aa90
+begin
+	x1lims_ackley = (-10, 10) #dimension 1
+	x2lims_ackley= (-10, 10)  #dimension 2
+end
+
+# ╔═╡ 680883a2-493f-11eb-0d58-03d0b5b302e2
+md"""The function can be run as easily as shown in the next block:"""
+
+# ╔═╡ 7974813e-493f-11eb-34ca-bf62397015eb
+begin
+	#initialize population
+	population = init_nests(15, x1lims_ackley, x2lims_ackley)
+
+	#run method
+	cuckoo!(ackley, population, x1lims_ackley, x2lims_ackley) 
+end
+
+# ╔═╡ f936c1f2-493f-11eb-2e6d-8bc27014712b
+md"""### Visualization"""
+
+# ╔═╡ dcd0b396-43af-11eb-1aa9-d3baa5317357
+md"""Using the implemented function `plot_solution` it is possible to visualize how by starting from the initial population (green points) the method reaches the final solution (white point):"""
+
+# ╔═╡ 279492ea-405b-11eb-2af7-13ec8531169a
+function plot_solution_2D(f, x1lims, x2lims)	
+	#initialize population
+	population = init_nests(20, x1lims, x2lims)
+
+	#draw function landscape
+	heatmap(x1lims[1]:0.1:x1lims[2], x2lims[1]:0.1:x2lims[2], f)
+
+
+	#plot initial population
+	for nest in population 
+    	scatter!([nest[1]], [nest[2]], color=:green, label="",
+            markersize=3)
+	end 
+
+	#run method
+	solution_stat = cuckoo!(f, population, x1lims, x2lims) 
+
+	#plot the final solution
+	s1=round(solution_stat[1][1],digits=1)	
+	s2=round(solution_stat[1][2],digits=1)
+	scatter!([solution_stat[1][1]],[solution_stat[1][2]], color=:white, 
+		     label="Solution [$(s1),$(s2)]", markersize=4)
+
+end
+
+# ╔═╡ c9c0363e-493e-11eb-25c5-271196c1a76f
+plot_solution_2D(ackley, x1lims_ackley, x2lims_ackley)
+
+# ╔═╡ 3559598a-447a-11eb-38cf-ad3f603b7063
+md"""We can also take a look at how the fitness improves at each round of optimization:"""
+
+# ╔═╡ fd697cfc-405b-11eb-38c8-91a9ad3ba93a
+
+
+# ╔═╡ 26deebee-446f-11eb-3390-4748fcdf5a36
+begin 
+	#initialize population
+	fitness_population = init_nests(15, x1lims_ackley, x2lims_ackley)  
+
+	fitness=Vector{Float64}(undef, 15)
+	#run the method for one generation 15 times to plot the fitness
+	for i=1:15
+		solution=cuckoo!(ackley, fitness_population, x1lims_ackley, x1lims_ackley, gen=1)
+		fitness[i]=solution[2]
+
+	end  
+	plot([fitness], label="", xlabel="Iteration Number", ylabel="Fitness")
+
+end
+
+# ╔═╡ 6a1c79aa-493d-11eb-09db-a10a40ca8094
+md"""### Animation"""
+
+# ╔═╡ 34675522-43b0-11eb-3d00-c12b0cf63781
+md"""In this animation we see how the solutions are updated at every round of optimization to reach the optimal solution, where the best solution is denoted by the white dot at each iteration."""
+
+# ╔═╡ 6ee961fc-407b-11eb-104c-4fd3d89b4790
+begin 	
+	#initialize population
+	moving_population = init_nests(15, x1lims_ackley, x2lims_ackley)  
+
+	#run the method for one generation 15 times to plot 15 population states
+	anim = @animate for i=1:15
+		solution=cuckoo!(ackley,moving_population,x1lims_ackley,x2lims_ackley, gen=1)
+
+		#plot landscape
+		heatmap(x1lims_ackley[1]:0.01:x1lims_ackley[2], 
+				x2lims_ackley[1]:0.01:x2lims_ackley[2], 
+				ackley)
+
+		for nest in moving_population
+			#population at iteration i
+			scatter!([nest[1]], [nest[2]], color=:green, label="",
+					markersize=3)
+		end
+		#best solution at iteration i
+		scatter!([solution[1][1]], [solution[1][2]], color=:white, label="", markersize=4)
+
+
+
+	end 
+
+gif(anim, fps = 1)
+end
+
+# ╔═╡ 43539fba-493d-11eb-0430-fbb7c588cba9
+md"""### Parameters and population size"""
+
+# ╔═╡ 857d667e-43b0-11eb-02c4-f97345274afd
+md"""The method allows to specify some parameters, namely:
+* `gen`: number of generation. If the number is too small you may get a suboptimal solution due to the fact that the problem may have not yet converged to its optimal solution. 
+* `alpha`: stepsize when transitioning to a new solution. By default it is 1.0 but can be varied depending on the dimension of the landscape.
+* `Pa`: rate of cuckoo's egg discovery. This parameter allows to balance how much the method relies on exploring known solutions (default `Pa`=0.25) as compared to inspecting new solution through Lévy flights.
+* `lambda`: exponent of the Lévy distribution used to sample stepsizes."""
+
+# ╔═╡ 703cb04e-446e-11eb-2954-ebf4cd33851e
+md"""Yet another element that strongly conditions the outcome of the method is the population size:"""
+
+# ╔═╡ 238232a4-4470-11eb-36b9-dfd43cc99640
+begin
+	plot()
+	populations = [2,5,10,20] 
+
+	for x = 1:length(populations) 
+		fitness=Vector{Float64}(undef, 15)			
+		fitness_population = init_nests(populations[x], x1lims_ackley, x2lims_ackley)  
+		#run the method for one generation 15 times to plot the fitness
+		for i=1:15
+			solution=cuckoo!(ackley, fitness_population, x1lims_ackley, x2lims_ackley, gen=1)	
+			fitness[i]= solution[2]
+		end  
+		plot!([fitness], label="Population size: $(populations[x])", xlabel="Number of iterations", ylabel="Fitness")
+	end
+	current()
+end
+
+# ╔═╡ 9f1bb16c-4941-11eb-2db4-87331f69ca7f
+md"""### Using more than 2 dimensions"""
+
+# ╔═╡ 6464ddae-4942-11eb-32ab-c55e665f41c6
+md"""The method can take a problem of any dimension, assumed that the accurate number of limits is provided. Here the 3-dimensional ackley function, having as global minimizer (0, 0, 0) is solved: """
+
+# ╔═╡ 24b5909a-4942-11eb-2768-a1586788bff2
+begin 
+	x1lims_3D = (-10, 10) #dimension 1
+	x2lims_3D= (-10, 10)  #dimension 2
+ 	x3lims_3D= (-10, 10)  #dimension 3
+
+	#initialize population
+	population_3D = init_nests(15, x1lims_3D, x2lims_3D, x3lims_3D)
+
+	#run method
+	cuckoo!(ackley, population_3D, x1lims_3D, x2lims_3D, x3lims_3D) 
+end
+
+# ╔═╡ 9a9b2b30-493d-11eb-2db1-f39a0782c896
+md"""### More complex functions"""
+
+# ╔═╡ cb3d1354-493b-11eb-29af-4d36cb31f4e4
+begin
+	rosenbrock((x1, x2); a=1, b=5) = (a-x1)^2 + b*(x2-x1^2)^2
+    rosenbrock(x1, x2; kwargs...) = rosenbrock((x1, x2); kwargs...)
+end
+
+# ╔═╡ 93c7588c-4eb1-11eb-1e5a-a7bedb46fefd
+md"""**Rosenbrock function** global minimum: (1,1)"""
+
+# ╔═╡ f1110a4e-493c-11eb-234d-ddc38917505c
+begin
+	x1lims_rastrigine = (-2.048, 2.048) #dimension 1
+	x2lims_rastrigine= (-2.048, 2.048)  #dimension 2
+	plot_solution_2D(rosenbrock, x1lims_rastrigine, x2lims_rastrigine)
+end
+
+# ╔═╡ a946be48-4eaf-11eb-1e20-a152a2278921
+begin
+	function branin((x1, x2); a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)
+    	return a * (x2 - b * x1^2 + c * x1 - r)^2 + s * (1 - t) * cos(x1) + s
+	end
+
+	branin(x1, x2; kwargs...) = branin((x1, x2); kwargs...)
+end
+
+# ╔═╡ 362f0eae-4eb1-11eb-26a2-15c63cfa9cf0
+md""" **Branin function** global minima: 
+* (-3.14, 12.275)
+* (3.14, 2.275)
+* (9.42,2.475)"""
+
+# ╔═╡ b64b72c2-4eb0-11eb-2434-1b0824d7fa17
+begin
+	x1lims_branin = (-5, 10) #dimension 1
+	x2lims_branin= (0, 15)  #dimension 2
+	plot_solution_2D(branin, x1lims_branin, x2lims_branin)
+end
+
+# ╔═╡ Cell order:
+# ╟─d03cf1ca-4059-11eb-3339-99c511eba3c8
+# ╟─db128a7a-4085-11eb-08b9-5b5199efecc3
+# ╟─3c8f220e-4086-11eb-175b-f92b3e3e8253
+# ╠═694cf0e0-405a-11eb-08f9-bd9a2e0c6232
+# ╟─56c0a760-4086-11eb-1c5f-45828eb40989
+# ╟─be3ef176-4086-11eb-365c-e1747053c12d
+# ╟─e9040d5a-4087-11eb-03a6-51e941460bc5
+# ╠═e0165c02-405a-11eb-296f-05e437e873e3
+# ╟─2607802e-4088-11eb-1a50-bbabd9cbeaee
+# ╠═6cfa745e-408c-11eb-345c-29b71c84aa90
+# ╟─680883a2-493f-11eb-0d58-03d0b5b302e2
+# ╠═7974813e-493f-11eb-34ca-bf62397015eb
+# ╟─f936c1f2-493f-11eb-2e6d-8bc27014712b
+# ╟─dcd0b396-43af-11eb-1aa9-d3baa5317357
+# ╠═c9c0363e-493e-11eb-25c5-271196c1a76f
+# ╟─279492ea-405b-11eb-2af7-13ec8531169a
+# ╟─3559598a-447a-11eb-38cf-ad3f603b7063
+# ╟─fd697cfc-405b-11eb-38c8-91a9ad3ba93a
+# ╟─26deebee-446f-11eb-3390-4748fcdf5a36
+# ╟─6a1c79aa-493d-11eb-09db-a10a40ca8094
+# ╟─34675522-43b0-11eb-3d00-c12b0cf63781
+# ╟─6ee961fc-407b-11eb-104c-4fd3d89b4790
+# ╟─43539fba-493d-11eb-0430-fbb7c588cba9
+# ╟─857d667e-43b0-11eb-02c4-f97345274afd
+# ╟─703cb04e-446e-11eb-2954-ebf4cd33851e
+# ╟─238232a4-4470-11eb-36b9-dfd43cc99640
+# ╟─9f1bb16c-4941-11eb-2db4-87331f69ca7f
+# ╟─6464ddae-4942-11eb-32ab-c55e665f41c6
+# ╠═24b5909a-4942-11eb-2768-a1586788bff2
+# ╟─9a9b2b30-493d-11eb-2db1-f39a0782c896
+# ╟─cb3d1354-493b-11eb-29af-4d36cb31f4e4
+# ╟─93c7588c-4eb1-11eb-1e5a-a7bedb46fefd
+# ╟─f1110a4e-493c-11eb-234d-ddc38917505c
+# ╟─a946be48-4eaf-11eb-1e20-a152a2278921
+# ╟─362f0eae-4eb1-11eb-26a2-15c63cfa9cf0
+# ╠═b64b72c2-4eb0-11eb-2434-1b0824d7fa17

project ideas.md

@@ -11,7 +11,7 @@ Optimization methods not covered in the theory. Provide a working implementation
 - Conjugate gradient descent
 - [Cross-entropy method](https://en.wikipedia.org/wiki/Cross-Entropy_Method)
 - Firefly algorithm
-- Bee Colony Optimization
+- ~~Bee Colony Optimization~~
 - Cuckoo search
 - maximum flow problems
 - Levenberg–Marquardt algorithm (for nonlinear regression and inverse kinematics)
@@ -23,7 +23,7 @@ Solve a problem using a method from the course or a new one.
 
 - time table scheduling
 - analytic solution of a ODE using genetic programming
-- solve a Sudoku using local search
+- ~~solve a Sudoku using local search~~
 - graph matching using simulated annealing
 - finding negative cycles using the Bellman–Ford algorithm
 - multiple sequence alignment using Genetic Algorithms

src/STMOZOO.jl

@@ -3,6 +3,8 @@ module STMOZOO
 
 # execute your source file and export the module you made
 include("example.jl")
-export Example
+include("cuckoo.jl")
+
+export Example, Cuckoo
 
 end # module