diff --git a/.gitignore b/.gitignore index cfb93320..7ad65f3d 100644 --- a/.gitignore +++ b/.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/ diff --git a/Project.toml b/Project.toml index 881f78bc..e3dc96b9 100644 --- a/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index 21d4327e..bfceeb44 100644 --- a/docs/make.jl +++ b/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" ]) #= diff --git a/docs/src/man/cuckoo.md b/docs/src/man/cuckoo.md new file mode 100644 index 00000000..6c91eee6 --- /dev/null +++ b/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! +``` \ No newline at end of file diff --git a/instructions.md b/instructions.md index 3fefa32e..00eb31d6 100644 --- a/instructions.md +++ b/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. diff --git a/notebook/cuckoo.jl b/notebook/cuckoo.jl new file mode 100644 index 00000000..8c8f1e97 --- /dev/null +++ b/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 diff --git a/project ideas.md b/project ideas.md index e0691d1a..ee99c8c2 100644 --- a/project ideas.md +++ b/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 diff --git a/src/STMOZOO.jl b/src/STMOZOO.jl index 34b94d12..5189a2a6 100644 --- a/src/STMOZOO.jl +++ b/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 diff --git a/src/cuckoo.jl b/src/cuckoo.jl new file mode 100644 index 00000000..48470fad --- /dev/null +++ b/src/cuckoo.jl @@ -0,0 +1,232 @@ +module Cuckoo + + # export functions relevant for the user + export init_nests, cuckoo! + + # import external packages + using SpecialFunctions + + """ + Nest + + Hold each solution and its fitness. + """ + mutable struct Nest + position::Array{Float64} + fitness::Float64 + end + + + """ + init_nests(n::Int, lims...) + + Create an array of `n` random solutions. + + Each dimension of the solution is constrained by an upper and a lower bound provided in + the `lims` parameter which automatically collects in a tuple all the couples of bounds + passed for each dimension. + + # Examples + + Create a population of 5 where each point has 2 dimensions, respectively having as lower + and upper bound `x1lims = (-10, 10)` and `x2lims=(-15, 15)`. + + ```julia + julia> x1lims=(-10,10) + julia> x2lims=(-15,15) + julia> p=init_nests(5,x1lims,x2lims) + 5-element Array{Array{Float64,1},1}: + [3.2506269365826235, -13.904103441708813] + [7.6966755817918155, 0.019170453403416943] + [9.177999793081582, -7.738154316325958] + [1.0837999423638358, 14.655548441822855] + [6.640425508361417, -4.823307495174905] + ``` + """ + function init_nests(n::Int, lims...) + @assert all([length(lim)==2 for lim in lims]) "Every bound should be composed of a lower and an upper limit." + @assert all([l <= u for (l, u) in lims]) "Every bound should be defined with the lower one as first element and the upper one as second" + @assert n>1 "The population should be composed of at least two individuals" + + return [([(u - l) * rand() + l for (l, u) in lims]) for i in 1:n] + end + + + """ + levy_step(d::Int, lambda::AbstractFloat) + + Return the step to be added to an existing solution. + + The step direction is randomly picked from a uniform distribution and the stepsize is + drawn from a Lévy distribution with exponent `lambda`, obtained through implementing a + simplified version of Mantegna algorithm. This allows to simulate a Lévy flight + performed by a cuckoo to reach a new nest. The parameter `d` indicates the number of + dimensions of the solution. + """ + function levy_step(d::Int, lambda::AbstractFloat) + #choice of random direction + dir = [rand() for x in 1:d] + + #generation of stepsize + sigma = (gamma(1 + lambda) * sin(pi * lambda/2)/ + (gamma((1 + lambda) / 2) * lambda * 2^((lambda - 1)/2)))^(1 / lambda) + + steps = Vector{Float64}(undef, d) + for x in 1:d + u = sigma*randn() + v = randn() + steps[x] = u ./ abs(v) .^ (1/lambda) + end + + return steps .* dir + end + + + """ + get_random_nest_index(n::Int; not::Set{Int}=Set([0])) + + Return a random index between 1 and `n`, representing one of the nests. + + If `not` is specified, the returned index must be different from any element in it. + """ + function get_random_nest_index(n::Int; not::Set{Int}=Set([0])) + @assert length(not) lims[d][2], Float64(lims[d][2]), new_pos[d])) for d in 1:n] + return new_pos + end + + + """ + cuckoo(f::Function, population::Array, lims...; + gen::Int=40, + Pa::AbstractFloat=0.25, + alpha::AbstractFloat=1.0, + lambda::AbstractFloat=1.5) + + Implement the backbone of the cuckoo search method and return a tuple containing the + solution and its fitness. + + At each iteration these two optimization steps are carried out: + + - Starting from the egg (solution) contained in each nest a Lévy flight is simulated to create a new solution. To preserve the total number of nests, this egg takes the place of a randomly chosen nest, but only if it has a better fitness. This step allows to diversify the exploration of the search space by performing a farfield optimization improved by the fact that the Lévy distribution is an heavy-tailed distribution; + + - A fraction `Pa` of worst eggs (bad solution) gets discovered by the host bird and the nest is abandoned. A new nest with a new solution is generated, biased toward two good quality eggs randomly picked from the population. This step allows to perform more of a locally intensified search by exploiting the neighborhood of current solutions. + + # Arguments + - `f:Function` function to minimize. + - `population::Array` array obtained from running the function [`init_nests`](@ref). + - `lims...` a number of parameters corresponding to the problem dimensionality, each is a couple (arrays/tuples) of lower and upper bounds. + - `gen::Int=40` number of generations. + - `Pa::AbstractFloat=0.25` rate of cuckoo eggs which gets discovered and abandoned at each generation. + - `alpha::AbstractFloat=1.0` a positive scaling parameter for step size. + - `lambda::AbstractFloat=1.5` the exponentiation parameter for the Lévy distribution. + + # Examples + ```julia + julia> x1lims = (-10,10) + julia> x2lims = (-15,15) + julia> population = init_nests(5,x1lims,x2lims) + julia> cuckoo!(function_name, population, x1lims, x2lims, gen=50, Pa=0.5, alpha=0.5, lambda=2.0) + ([-0.00401120574383079, 0.00046356024376423543], -22.706426805416967) + ``` + """ + function cuckoo!(f::Function, population::Array, lims...; + gen::Int=40, + Pa::AbstractFloat=0.25, + alpha::AbstractFloat=1.0, + lambda::AbstractFloat=1.5) + @assert (lambda<=3.0 && lambda>=1.0) "Lambda should be between 1 and 3 (included)" + @assert (alpha>0.0) "Alpha should be a positive value" + @assert (Pa<=1.0 && Pa>=0.0) "Pa should be a value between 0 and 1 (included)" + + n = length(population) #number of nests + d = length(lims) #problem dimensionality + + #create array of Nest structs to hold both the position and the fitness + nests = Nest.(population, f.(population)) + + #save current best solution + best_index = argmin([nests[i].fitness for i in 1:n]) + best_pos = nests[best_index].position + best_fit = nests[best_index].fitness + + #main loop for number of generations + while gen > 0 + #from each nest a cuckoo performs a flight to a new position + for i in 1:n + #generate a new solution with Lévy flights + new_pos = nests[i].position .+ alpha*levy_step(d, lambda) + new_pos = check_limits(new_pos, lims...) + new_fit = f(new_pos) + + #choose a nest randomly, different from i + j = get_random_nest_index(n, not=Set(i)) + old_fit = nests[j].fitness + + #replace j by new solution i if the fitness improves + if new_fit < old_fit + nests[j].position = new_pos + nests[j].fitness = new_fit + end + + end + + #sort from worst to best solution + sort!(nests, by = n -> n.fitness, rev=true) + + #a fraction Pa of worst nests (high fitness) gets discovered by the host bird + #and abandoned (at least two nests are not discovered) + n_worst = floor(Int, n * Pa) + n_worst = n_worst > n-2 ? n-2 : n_worst + + i = 1 + while i < n_worst + #randomly pick two (different) good solutions + j = get_random_nest_index(n, not=Set(1:n_worst)) + k = get_random_nest_index(n, not=push!(Set(1:n_worst),j)) + nj = nests[j].position + nk = nests[k].position + + #new solution generated by combining the two random good solutions + new_pos = nests[i].position .+ alpha.*rand().*(nj .- nk) + new_pos = check_limits(new_pos, lims...) + new_fit = f(new_pos) + + nests[i].position = new_pos + nests[i].fitness = new_fit + + i +=1 + end + + #save best solution + best_index = argmin([nests[i].fitness for i in 1:n]) + best_pos = nests[best_index].position + best_fit = nests[best_index].fitness + + #for visualization purposes + population .= [nests[i].position for i in 1:n] + + gen -= 1 + end + return best_pos, best_fit + end +end \ No newline at end of file diff --git a/test/cuckoo.jl b/test/cuckoo.jl new file mode 100644 index 00000000..10ca8bca --- /dev/null +++ b/test/cuckoo.jl @@ -0,0 +1,102 @@ +# Unit test on the Cuckoo module + +@testset "Cuckoo" begin + using STMOZOO.Cuckoo + using STMOZOO.Cuckoo: levy_step, get_random_nest_index, check_limits + + @testset "initialize nests" begin + n = 5 + xlims=((-10,10), (-20,20), (-10,10)) + extr_lims=((-10,-10),(-10,-10)) + wrong_lims_1=((10,-10), (-10,10)) + wrong_lims_2=(5, (-10,10)) + wrong_lims_3=((-5,5,5), (-10,10)) + + #type test + @test init_nests(n, xlims[1], xlims[2], xlims[3], xlims[3]) isa Array + @test init_nests(n, extr_lims[1], extr_lims[2]) isa Array + @test init_nests(2, xlims[1], xlims[2], xlims[3]) isa Array + + #dimension test + @test length(init_nests(n, xlims[1], xlims[2], xlims[3])) == n + + #correctness first point + @test all([point[1] >= xlims[1][1] for point in init_nests(n, xlims[1], xlims[2], xlims[3])]) + @test all([point[1] <= xlims[1][2] for point in init_nests(n, xlims[1], xlims[2], xlims[3])]) + + #correctness last point + @test all([point[length(xlims)] >= last(xlims)[1] for point in init_nests(n, xlims[1], xlims[2], xlims[3])]) + @test all([point[length(xlims)] <= last(xlims)[2] for point in init_nests(n, xlims[1], xlims[2], xlims[3])]) + + #assertion errors + @test_throws AssertionError init_nests(n, wrong_lims_1[1], wrong_lims_1[2]) + @test_throws AssertionError init_nests(n, wrong_lims_2[1], wrong_lims_2[2]) + @test_throws AssertionError init_nests(n, wrong_lims_3[1], wrong_lims_3[2]) + @test_throws AssertionError init_nests(1, xlims[1], xlims[2], xlims[3]) + end + + d = 3 + lambda= 3/2 + @test levy_step(d, lambda) isa Array + + @testset "Get random nest index" begin + n=3 + + #type test + @test get_random_nest_index(n) isa Int + @test get_random_nest_index(n, not=Set([1])) isa Int + + #correctness + @test get_random_nest_index(n, not=Set([1,2])) == 3 + + #assertion errors + @test_throws AssertionError get_random_nest_index(n, not=Set([1,2,3])) + + end + + @testset "Check limits" begin + wrong_pos = [1.23, 2.34, -1.2] + pos_in=[5.0, 5.0] + pos_out_down_one=[-11.0, 5.0] + pos_out_up_two=[5.0, 11.0] + pos_out_all=[-11.0,11.0] + xlims=((-10,10),(-10,10)) + + #type test + @test check_limits(pos_in, xlims[1], xlims[2]) isa Array{Float64} + + #correctness + @test check_limits(pos_in, xlims[1], xlims[2]) == pos_in + @test check_limits(pos_out_down_one, xlims[1], xlims[2]) == [-10.0,5.0] + @test check_limits(pos_out_up_two, xlims[1], xlims[2]) == [5.0,10.0] + @test check_limits(pos_out_all, xlims[1], xlims[2]) == [-10.0,10.0] + + #assertion errors + @test_throws AssertionError check_limits(wrong_pos, xlims[1], xlims[2]) + end + + @testset "Cuckoo method" begin + #create function for testing the method + 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 + + #initialize population + x1lims = (-10, 10) + x2lims = (-10, 10) + population = init_nests(25, x1lims, x2lims) + + #type tests + @test cuckoo!(ackley, population, x1lims, x2lims) isa Tuple + @test cuckoo!(ackley, population, x1lims, x2lims)[1] isa Array + @test cuckoo!(ackley, population, x1lims, x2lims)[2] isa AbstractFloat + + #assertion errors + @test_throws AssertionError cuckoo!(ackley, population, x1lims, x2lims, lambda=4.0) + @test_throws AssertionError cuckoo!(ackley, population, x1lims, x2lims, alpha=-1.0) + @test_throws AssertionError cuckoo!(ackley, population, x1lims, x2lims, Pa=4.0) + + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index ff9b0f7e..8fd3222c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ using Test include("example.jl") -# add here the file with your unit tests \ No newline at end of file +include("cuckoo.jl") \ No newline at end of file