diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 16f796cc..125cd3ef 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -4,7 +4,8 @@ on: branches: - 'master' pull_request: - types: [opened, reopened] + +jobs: jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} @@ -13,7 +14,7 @@ jobs: fail-fast: false matrix: version: - - 1.5 + - 1.7 os: - ubuntu-latest arch: diff --git a/.github/workflows/documentation-github-actions.yml b/.github/workflows/documentation-github-actions.yml index 828d10e0..378ebba7 100644 --- a/.github/workflows/documentation-github-actions.yml +++ b/.github/workflows/documentation-github-actions.yml @@ -14,7 +14,7 @@ jobs: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@latest with: - version: '1.5' + version: '1.7' - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy diff --git a/Project.toml b/Project.toml index 67e9c7e5..c5d1051a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Metaheuristics" uuid = "bcdb8e00-2c21-11e9-3065-2b553b22f898" authors = ["Jesus Mejia "] -version = "3.2.16" +version = "3.3.0" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" @@ -10,16 +10,20 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +SearchSpaces = "eb7571c6-2196-4f03-99b8-52a5a35b3163" SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] Distances = "^0.10.7" JMcDM = "^0.7.1" +Reexport = "^1" Requires = "1" +SearchSpaces = "^0.2" SnoopPrecompile = "1" -julia = "^1.5" +julia = "^1.7" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" diff --git a/README.md b/README.md index 2d8f319e..051219af 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ High-performance metaheuristics for global optimization. ## Installation -Open the Julia REPL (v1.1 or later) and press `]` to open the Pkg prompt. To add this package, use the add command: +Open the Julia REPL and press `]` to open the Pkg prompt. To add this package, use the add command: ``` pkg> add Metaheuristics @@ -105,15 +105,16 @@ Code the objective function: f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x) ) ``` -Instantiate the bounds, note that `bounds` should be a $2\times 10$ `Matrix` where -the first row corresponds to the lower bounds whilst the second row corresponds to the -upper bounds. +Instantiate the bounds: ```julia D = 10 -bounds = [-5ones(D) 5ones(D)]' +bounds = boxconstraints(lb = -5ones(D), ub = 5ones(D)) ``` +Also, `bounds` can be a $2\times 10$ `Matrix` where the first row corresponds to the +lower bounds whilst the second row corresponds to the upper bounds. + Approximate the optimum using the function `optimize`. ```julia diff --git a/docs/Project.toml b/docs/Project.toml index 4bea2e4a..c7d98d08 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Documenter = "0.27.15" diff --git a/docs/make.jl b/docs/make.jl index 99a1a041..07f66f7f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,7 +9,8 @@ makedocs( prettyurls = get(ENV, "CI", nothing) == "true", assets = ["assets/favicon.ico", "assets/extra_styles.css"], analytics = "UA-184071594-1", - collapselevel = 2, + collapselevel = 1, + ansicolor=true, ), sitename="Metaheuristics.jl", authors = "Jesús Mejía", @@ -19,6 +20,7 @@ makedocs( "tutorials/simple-tutorial.md", "tutorials/create-metaheuristic.md", "tutorials/parallelization.md", + "tutorials/n-queens.md", ], "Examples" => "examples.md", "Algorithms" => "algorithms.md", @@ -27,6 +29,7 @@ makedocs( "Multi-Criteria Decision Making" => "mcdm.md", "Visualization" => "visualization.md", "API References" => "api.md", + "FAQ" => "faq.md", "Contributing" => "contributing.md", "References" => "references.md", ] diff --git a/docs/src/api.md b/docs/src/api.md index 5a0ef53b..b5e3ed7d 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -198,6 +198,10 @@ GenerationalReplacement ElitistReplacement ``` +```@docs +RankAndCrowding +``` + ## Population diff --git a/docs/src/assets/extra_styles.css b/docs/src/assets/extra_styles.css index 381c4033..48e258fe 100644 --- a/docs/src/assets/extra_styles.css +++ b/docs/src/assets/extra_styles.css @@ -4,3 +4,39 @@ dl dd div:target { border-radius: 10px; } +html.theme--documenter-dark #documenter .docs-main, +#documenter .docs-main{ + max-width: 45rem; +} + + +pre { + border-radius: 7px; +} + +article.docstring { + border:none; + box-shadow: 1px 1px 9px -2px rgba(0,0,0,0.7); + -webkit-box-shadow: 1px 1px 9px -2px rgba(0,0,0,0.7); + -moz-box-shadow: 1px 1px 9px -2px rgba(0,0,0,0.7); + margin-bottom: 2em; +} + +code.nohighlight.ansi { + transition: opacity .5s ease-out; + opacity: 0.5; +} + +.documenter-example-output code.nohighlight.ansi{ + opacity: 0.8; +} + + +pre:hover > code.nohighlight.ansi { + opacity: 9; +} + +.content pre { + border: 1px solid #e6e6e6; +} + diff --git a/docs/src/examples.md b/docs/src/examples.md index 33d92e99..9b321769 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -6,29 +6,66 @@ examples for using the implemented `Metaheuristics`. ## Single-Objective Optimization -```@repl +Firstly import this package + +```@example SingleObjective using Metaheuristics +``` +Now, let us define the objective function to be minimized: -f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x) ) # objective function +```@example SingleObjective +f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x) ) +``` -bounds = [-5ones(10) 5ones(10)]' # limits/bounds +The search space (a.k.a. box-constraints) can be defined as follows: -information = Information(f_optimum = 0.0); # information on the minimization problem +```@example SingleObjective +bounds = boxconstraints(lb = -5ones(10), ub = 5ones(10)) +``` -options = Options(f_calls_limit = 9000*10, f_tol = 1e-5); # generic settings +!!! compat "boxconstraints in a Matrix format." + You can also define the bounds using `bounds = [-5ones(10) 5ones(10)]'`; however this + is not longer recommended. -algorithm = ECA(information = information, options = options) # metaheuristic used to optimize +It is possible to provide some information on the minimization problem. +Let's provide the true optimum to stop the optimizer when a tolerance `f_tol` is satisfied. -result = optimize(f, bounds, algorithm) # start the minimization process +```@example SingleObjective +information = Information(f_optimum = 0.0) +``` +Generic options or settings (e.g. budget limitation, tolerances, etc) can be provided as follows: + +```@example SingleObjective +options = Options(f_calls_limit = 9000*10, f_tol = 1e-5, seed=1) +``` + +Now, we can provide the Information and Options to the optimizer (ECA in this example). + +```@example SingleObjective +algorithm = ECA(information = information, options = options) +``` +Now, the optimization is performed as follows: + +```@example SingleObjective +result = optimize(f, bounds, algorithm) +``` + +The minimum and minimizer: + +```@example SingleObjective minimum(result) +``` + +```@example SingleObjective minimizer(result) +``` -result = optimize(f, bounds, algorithm) # note that second run is faster +!!! compat "Second run is faster" + As you may know, the second run can be faster: -``` ## Constrained Optimization @@ -81,6 +118,7 @@ you need to provide constraints if they exist, otherwise put `gx = [0.0]; hx = [ to indicate an unconstrained multiobjective problem. ```@repl +using UnicodePlots # to visualize in console (optional) using Metaheuristics function f(x) @@ -244,7 +282,7 @@ algorithm performance or test new mechanisms. This example illustrates how to do Let's assume that we want to modify the stop criteria for `ECA`. See [Contributing](@ref) for more details. -```@repl +```julia using Metaheuristics import LinearAlgebra: norm diff --git a/docs/src/faq.md b/docs/src/faq.md new file mode 100644 index 00000000..eb01b40b --- /dev/null +++ b/docs/src/faq.md @@ -0,0 +1,10 @@ +# FAQ + +Answers to Frequently Asked Questions: + +## How to solve combinatorial problems? + +This package was initially developed for numerical optimization, but recent updates +can handle combinatorial problems. The Genetic Algorithm framework ([`GA`](@ref)) can be used. + + diff --git a/docs/src/index.md b/docs/src/index.md index 8fb2c2a3..d2f56d0f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,6 @@ # Metaheuristics - an Intuitive Package for Global Optimization -**Author: Jesus Mejía (@jmejia8)** +**Author: Jesús-Adolfo Mejía-de-Dios (@jmejia8)** High-performance algorithms for optimization coded purely in a high-performance language. @@ -26,13 +26,13 @@ global optimization. Open the Julia (Julia 1.1 or Later) REPL and press `]` to open the Pkg prompt. To add this package, use the add command: -``` +```julia-repl pkg> add Metaheuristics ``` Or, equivalently, via the `Pkg` API: -```julia +```julia-repl julia> import Pkg; Pkg.add("Metaheuristics") ``` @@ -53,36 +53,45 @@ dimension number, assume $D=10$. Firstly, import the Metaheuristics package: -```julia +```@example julia using Metaheuristics ``` Code the objective function: -```julia + +```@example julia f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x) ) +nothing # hide ``` -Instantiate the bounds, note that `bounds` should be a $2\times 10$ `Matrix` where -the first row corresponds to the lower bounds whilst the second row corresponds to the -upper bounds. +Instantiate the bounds: -```julia +```@example julia D = 10 -bounds = [-5ones(D) 5ones(D)]' +bounds = boxconstraints(lb = -5ones(D), ub = 5ones(D)) +nothing # hide ``` +Also, `bounds` can be a $2\times 10$ `Matrix` where the first row corresponds to the +lower bounds whilst the second row corresponds to the upper bounds. + Approximate the optimum using the function `optimize`. -```julia +```@example julia +import Random: seed! # hide +seed!(50) # hide result = optimize(f, bounds) ``` Optimize returns a `State` datatype which contains some information about the approximation. For instance, you may use mainly two functions to obtain such an approximation. -```julia -@show minimum(result) -@show minimizer(result) +```@example julia +minimum(result) +``` + +```@example julia +minimizer(result) ``` diff --git a/docs/src/tutorials/n-queens.md b/docs/src/tutorials/n-queens.md new file mode 100644 index 00000000..ceef3f7d --- /dev/null +++ b/docs/src/tutorials/n-queens.md @@ -0,0 +1,99 @@ +# N-Queens + +In the N-queen problem, $N$ queens must be placed on an $N\times N$ chessboard without interfering with each other. +There can never be more than one queen on a chess row, column or diagonal because of the queen's ability to move vertically, horizontally and diagonally. + +## Solution Representation + +A naive representation is to use a binary matrix to represent a chessboard $c_{ij}$ +where $c_{ij}=1$ would indicate a queen in row $i$ and column $j$; however the search +space cardinality (number of solutions) is huge for this representation. +A permutation-based representation can be used instead. +A permutation will contain in position $j$ the row $i$ where a queen is, removing a lot of +infeasible solutions (queens attacking horizontally and vertically). + +Let's compute the number of solutions (a.k.a. cardinality): + +```@repl +using Metaheuristics +N = 8 # for 8-queens +cardinality(BitArraySpace(N*N)) +cardinality(PermutationSpace(N)) +``` + +## Objective Function + +It is necessary to minimize the number of attacks, the following objective function is +proposed to this end: + +```@example queens +function attacks(chessboard) + N = length(chessboard) + n_attacks = 0 + # check attack in both diagonas for each queen + for i = 1:N + Δrows = i + chessboard[i] + Δcols = i - chessboard[i] + for j = (i+1):N + # check diagonal [\] + n_attacks += j + chessboard[j] == Δrows ? 1 : 0 + # check inverse diagonal [/] + n_attacks += j - chessboard[j] == Δcols ? 1 : 0 + end + end + 2n_attacks +end +``` + +It can observed that horizontal and vertical are not considered due to the +adopted solution representation. + + +## Let's find a solution + +To minimize the number of attacks, let's use a Genetic Algorithm: + +```@example queens +using Metaheuristics # hide +N = 8 +optimize(attacks, PermutationSpace(N), GA); # hide +result = optimize(attacks, PermutationSpace(N), GA) +``` + +It can be observed that the number of attacks is: + +```@example queens +n_attacks = minimum(result) +``` + +The optimal permutation is: + +```@example queens +perm = minimizer(result) +``` + +Corresponding chessboard configuration: + +```@example queens +chessboard = zeros(Bool, N, N) +for (i,j) = enumerate(perm); chessboard[i,j] = true;end +chessboard +``` + +## The 20-queens case + +```@example queens +N = 20 +perm = optimize(attacks, PermutationSpace(N), GA, seed=1) |> minimizer +``` + +```@example queens +attacks(perm) +``` + + +```@example queens +chessboard = zeros(Bool, N, N) +for (i,j) = enumerate(perm); chessboard[i,j] = true;end +chessboard +``` diff --git a/docs/src/tutorials/simple-tutorial.md b/docs/src/tutorials/simple-tutorial.md index 0317a65f..c71342bf 100644 --- a/docs/src/tutorials/simple-tutorial.md +++ b/docs/src/tutorials/simple-tutorial.md @@ -16,21 +16,22 @@ Note that the global optimum is obtained when $x_i = 0$ for all $i$. Thus, $\min Objective function: -```julia +```@example julia +using Metaheuristics # hide f(x) = 10length(x) + sum( x.^2 - 10cos.(2pi*x) ) ``` Bounds: -```julia -bounds = [-5ones(10) 5ones(10)]' +```@example julia +bounds = BoxConstrainedSpace(lb = -5ones(10), ub = 5ones(10)) ``` ## Providing Information Since the optimum is known, then we can provide this information to the optimizer. -```julia +```@example julia information = Information(f_optimum = 0.0) ``` @@ -41,8 +42,9 @@ of function evaluations. To do that, let's assume that the metaheuristic should $9000D$ times the objective function. Moreover, since `information` is provided, then we can set the desired accuracy ($|f(x) - f(x^*)| $) to $10^{-5}$. -```julia -options = Options(f_calls_limit = 9000*10, f_tol = 1e-5) +```@example julia +options = Options(f_calls_limit = 9000*10, f_tol = 1e-5); +nothing # hide ``` ## Choose a Metaheuristic @@ -55,7 +57,7 @@ algorithm following the same steps. The metaheuristics accept their parameters but share two common and **optional** settings `information` and `options`. -```julia +```@example julia algorithm = ECA(information = information, options = options) ``` @@ -67,23 +69,23 @@ algorithm = ECA(information = information, options = options) Now, we are able to approximate the optimum. To do that it is necessary to use the `optimize` function as follows: -```julia +```@example julia result = optimize(f, bounds, algorithm) ``` ## Get the Results -Once `optimize` stops, we can get the approximate solutions. +Once [`optimize`](@ref) stops, we can get the approximate solutions. Approximated minimum: -```julia +```@example julia fx = minimum(result) ``` Approximated minimizer: -```julia +```@example julia x = minimizer(result) ``` @@ -97,7 +99,7 @@ get their positions. We recommend you wrap your program in a function for performance purposes: -```julia +```@example julia using Metaheuristics function main() @@ -105,7 +107,7 @@ function main() f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x) ) # limits/bounds - bounds = [-5ones(10) 5ones(10)]' + bounds = BoxConstrainedSpace(lb = -5ones(10), ub = 5ones(10)) # information on the minimization problem information = Information(f_optimum = 0.0) @@ -119,14 +121,13 @@ function main() # start the minimization process result = optimize(f, bounds, algorithm) - fx = minimum(result) x = minimizer(result) - @show fx - @show x + x, fx end +main() ``` ## Summary diff --git a/src/Metaheuristics.jl b/src/Metaheuristics.jl index ffdf5443..7b34ec3e 100644 --- a/src/Metaheuristics.jl +++ b/src/Metaheuristics.jl @@ -22,17 +22,17 @@ export GA, RandomInBounds, RandomBinary, RandomPermutation export TournamentSelection, RouletteWheelSelection export UniformCrossover, OrderCrossover, SBX export BitFlipMutation, SlightMutation, PolynomialMutation -export GenerationalReplacement, ElitistReplacement +export GenerationalReplacement, ElitistReplacement, RankAndCrowding export sample, LatinHypercubeSampling, Grid export εDE, Restart export optimize! export set_user_solutions! +export boxconstraints include("externals.jl") include("core/abstracts.jl") -include("core/stop_status_codes.jl") include("core/information.jl") include("core/options.jl") include("core/state.jl") @@ -45,11 +45,12 @@ include("solutions/display.jl") include("operators/operators.jl") +include("termination/termination.jl") + include("common/utils.jl") include("common/multi-objective-functions.jl") include("common/repair.jl") -include("common/stop.jl") include("common/compare.jl") include("common/non-dominated-sorting.jl") include("common/set_user_solutions.jl") @@ -61,8 +62,9 @@ include("PerformanceIndicators/PerformanceIndicators.jl") include("algorithms/algorithm.jl") +include("parameters/parameters.jl") -include("optimize.jl") +include("optimize/optimize.jl") ####################################################### # # diff --git a/src/PerformanceIndicators/covering.jl b/src/PerformanceIndicators/covering.jl index e9394cc1..e41ad06f 100644 --- a/src/PerformanceIndicators/covering.jl +++ b/src/PerformanceIndicators/covering.jl @@ -32,16 +32,19 @@ end covering(A::Vector{T}, B::Vector{T}) where T <: AbstractMultiObjectiveSolution = covering(fval.(A), fval.(B)) -function covering(A::State{T}, B::State{T}) where T <: AbstractMultiObjectiveSolution +function covering(A::State{T}, B::State{T}; + verbose = true + ) where T <: AbstractMultiObjectiveSolution + A_non_dominated = get_non_dominated_solutions(A.population) B_non_dominated = get_non_dominated_solutions(B.population) length(A_non_dominated) != length(A.population) && - @warn "Some solutions in B dominate each other (solutions ignored)." + verbose && @warn "Some solutions in B dominate each other (solutions ignored)." length(B_non_dominated) != length(B.population) && - @warn "Some solutions in B dominate each other (solutions ignored)." + verbose && @warn "Some solutions in B dominate each other (solutions ignored)." covering(A_non_dominated, B_non_dominated) end diff --git a/src/algorithms/ABC/ABC.jl b/src/algorithms/ABC/ABC.jl index 4f7f4218..9ab33f41 100644 --- a/src/algorithms/ABC/ABC.jl +++ b/src/algorithms/ABC/ABC.jl @@ -54,16 +54,14 @@ function ABC(; Ne = div(N+1, 2), No = div(N+1, 2), limit=10, - information = Information(), - options = Options() + kargs... ) parameters = ABC(N, Ne, No, limit) Algorithm( - parameters, - information = information, - options = options, + parameters; + kargs... ) @@ -82,7 +80,7 @@ function initialize!( options.parallel_evaluation && error("ABC is not supporting parallel evaluation. Put `options.parallel_evaluation=false`") - D = size(problem.bounds, 2) + D = getdim(problem) if options.f_calls_limit == 0 options.f_calls_limit = 10000D @@ -113,19 +111,19 @@ function update_state!( kargs... ) - D = size(problem.bounds, 2) + D = getdim(problem) fobj = problem.f bees = status.population Ne = parameters.Ne No = parameters.No - bounds = problem.bounds - a = view(bounds, 1,:) - b = view(bounds, 2,:) + bounds = problem.search_space + a = bounds.lb + b = bounds.ub - employedPhase!(bees,problem, Ne) - outlookerPhase!(bees,problem, No) + employedPhase!(bees,problem, Ne, options.rng) + outlookerPhase!(bees,problem, No, options.rng) - @inline genPos(D=D, a=Array(a), b = Array(b)) = a + (b - a) .* rand(D) + @inline genPos(D=D, a=a, b=b) = a + (b - a) .* rand(options.rng, D) best = chooseBest(bees, status.best_sol) status.f_calls += Ne + No + scoutPhase!(bees, problem, genPos, parameters.limit) diff --git a/src/algorithms/ABC/bee_dynamics.jl b/src/algorithms/ABC/bee_dynamics.jl index 6489408b..affded73 100644 --- a/src/algorithms/ABC/bee_dynamics.jl +++ b/src/algorithms/ABC/bee_dynamics.jl @@ -23,6 +23,8 @@ get_position(bee::Bee) = bee.sol.x Get the fitness of a bee when optimize using ABC algorithm. """ fval(bee::Bee) = bee.sol.f +fvals(bees::Vector{<:Bee}) = fval.(bees) +is_feasible(bee::Bee) = is_feasible(bee.sol) minimum(st::State{Bee{xf_indiv}}) = st.best_sol.sol.f @@ -48,14 +50,14 @@ function fit(fx) 1.0 + abs(fx) end -function updateBee!(bee, bee2, problem) +function updateBee!(bee, bee2, problem, rng) D = length(bee.sol.x) - ϕ = -1.0 + 2.0rand() + ϕ = -1.0 + 2.0rand(rng) v = ϕ*(bee.sol.x - bee2.sol.x) x_new = bee.sol.x + v - replace_with_random_in_bounds!(x_new, problem.bounds) + replace_with_random_in_bounds!(x_new, problem.search_space, rng) new_sol = create_solution(x_new, problem) @@ -68,18 +70,18 @@ function updateBee!(bee, bee2, problem) end end -function getK(i, N) - j = rand(1:N) +function getK(i, N, rng) + j = rand(rng, 1:N) while j == i - j = rand(1:N) + j = rand(rng, 1:N) end j end -function employedPhase!(bees, problem, Ne) +function employedPhase!(bees, problem, Ne, rng) N = length(bees) - for i in randperm(N)[1:Ne] - updateBee!(bees[i], bees[getK(i, N)], problem) + for i in randperm(rng, N)[1:Ne] + updateBee!(bees[i], bees[getK(i, N, rng)], problem, rng) end end @@ -97,13 +99,13 @@ function roulettSelect(bees, sum_f) return length(bees) end -function outlookerPhase!(bees, problem, No::Int) +function outlookerPhase!(bees, problem, No::Int, rng) N = length(bees) sum_f = sum(map(x->x.fit, bees)) for i=1:No j = roulettSelect(bees, sum_f) - updateBee!(bees[j], bees[getK(j, N)], problem) + updateBee!(bees[j], bees[getK(j, N, rng)], problem, rng) end end diff --git a/src/algorithms/CCMO/CCMO.jl b/src/algorithms/CCMO/CCMO.jl index a620af68..55711b72 100644 --- a/src/algorithms/CCMO/CCMO.jl +++ b/src/algorithms/CCMO/CCMO.jl @@ -74,13 +74,12 @@ stop reason: Maximum number of iterations exceeded. """ function CCMO(base_algorithm; - information = Information(), - options = Options(), + kargs... ) parameters = CCMO(base_algorithm.parameters, xFgh_indiv[],Float64[],Float64[]) - Algorithm(parameters, information = information, options = options) + Algorithm(parameters; kargs...) end function initialize!( @@ -150,10 +149,10 @@ function reproduction(status, parameters::CCMO{NSGA2}, problem) # selection _s = TournamentSelection(;N=base.N ÷ 2) # crossover - _c = SBX(η=base.η_cr, p = base.p_cr, bounds = problem.bounds) + _c = SBX(η=base.η_cr, p = base.p_cr, bounds = problem.search_space) # mutation - base.p_m < 0.0 && (base.p_m = 1/size(bounds,2)) - _m = PolynomialMutation(η=base.η_m, p = base.p_m, bounds = problem.bounds) + base.p_m < 0.0 && (base.p_m = 1/getdim(problem)) + _m = PolynomialMutation(η=base.η_m, p = base.p_m, bounds = problem.search_space) # population 1 parent_mask = selection(population, _s, fitness=fitness1) diff --git a/src/algorithms/CGSA/CGSA.jl b/src/algorithms/CGSA/CGSA.jl index 10a911e1..e73177de 100644 --- a/src/algorithms/CGSA/CGSA.jl +++ b/src/algorithms/CGSA/CGSA.jl @@ -99,17 +99,15 @@ function CGSA(; X::Matrix{Float64} = zeros(0,0), V::Matrix{Float64} = zeros(0,0), fitness::Vector{Float64} = zeros(0), - information = Information(), - options = Options() + kargs... ) parameters = CGSA(N, chValueInitial, chaosIndex, ElitistCheck, Rpower, Rnorm, wMax, wMin, X, V, fitness) Algorithm( - parameters, - information = information, - options = options, + parameters; + kargs... ) end @@ -126,12 +124,13 @@ function initialize!( Rnorm = parameters.Rnorm N = parameters.N - D = size(problem.bounds, 2) + D = getdim(problem) fobj = problem.f # bounds vectors - low, up = problem.bounds[1,:], problem.bounds[2,:] + low = problem.search_space.lb + up = problem.search_space.ub max_it = 500 options.iterations = options.iterations == 0 ? max_it : options.iterations @@ -199,17 +198,17 @@ function update_state!( #Calculation of accelaration in gravitational field. eq.7-10,21. - a = Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it) + a = Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it, options.rng) #Agent movement. eq.11-12 - X, V = move(X,a,V) + X, V = move(X,a,V, options.rng) # # Checking allowable range. # X = correctPop(X, low, up) for i = 1:N - x = reset_to_violated_bounds!(X[i,:], problem.bounds) + x = reset_to_violated_bounds!(X[i,:], problem.search_space) X[i,:] = x end diff --git a/src/algorithms/CGSA/physics.jl b/src/algorithms/CGSA/physics.jl index 20afcbe2..ef488f79 100644 --- a/src/algorithms/CGSA/physics.jl +++ b/src/algorithms/CGSA/physics.jl @@ -47,7 +47,7 @@ function Gconstant(iteration,max_it) end #This function calculates the accelaration of each agent in gravitational field. eq.7-10,21. -function Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it) +function Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it, rng) N, D = size(X) final_per= 2 #In the last iteration, only 2 percent of agents apply force to the others. @@ -68,7 +68,7 @@ function Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it) if j != i R = norm(X[i,:]-X[j,:], Rnorm) #Euclidian distanse. for k=1:D - E[i, k]= E[i,k]+rand() * (M[j]) * ((X[j,k]-X[i,k])/(R^Rpower + eps())) + E[i, k]= E[i,k]+rand(rng) * (M[j]) * ((X[j,k]-X[i,k])/(R^Rpower + eps())) #note that Mp(i)/Mi(i)=1 end end @@ -82,10 +82,10 @@ function Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it) end # This function updates the velocity and position of agents. -function move(X,a,V) +function move(X,a,V, rng) # movement. N, D = size(X) - V = rand(N, D).*V + a # eq. 11. + V = rand(rng, N, D).*V + a # eq. 11. X = X + V # eq. 12. return X, V diff --git a/src/algorithms/DE/DE.jl b/src/algorithms/DE/DE.jl index 64346ff7..0a24b67b 100644 --- a/src/algorithms/DE/DE.jl +++ b/src/algorithms/DE/DE.jl @@ -19,7 +19,7 @@ include("epsilonDE.jl") DE(; N = 0, F = 1.0, - CR = 0.9, + CR = 0.5, strategy = :rand1, information = Information(), options = Options() @@ -57,26 +57,21 @@ julia> optimize(f, [-1 -1 -1; 1 1 1.0], DE(N=50, F=1.5, CR=0.8)) """ function DE(; N::Int = 0, - F = 1.0, - CR = 0.9, + F = 0.7, + CR = 0.5, CR_min = CR, CR_max = CR, F_min = F, F_max = F, strategy::Symbol = :rand1, - information = Information(), - options = Options(), + kargs... ) parameters = DE(N, promote(F, CR, CR_min, CR_max, F_min, F_max)..., strategy) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm(parameters; kargs...) end @@ -96,14 +91,15 @@ function update_state!( CR = parameters.CR + rng = options.rng # stepsize if parameters.F_min < parameters.F_max - F = parameters.F_min + (F_max - parameters.F_min) * rand() + F = parameters.F_min + (F_max - parameters.F_min) * rand(rng) end if parameters.CR_min < parameters.CR_max CR = - parameters.CR_min + (parameters.CR_max - parameters.CR_min) * rand() + parameters.CR_min + (parameters.CR_max - parameters.CR_min) * rand(rng) end new_vectors = reproduction(status, parameters, problem) @@ -155,7 +151,7 @@ function initialize!( args...; kargs... ) - D = size(problem.bounds, 2) + D = getdim(problem) @@ -207,14 +203,14 @@ function reproduction(status, parameters::AbstractDifferentialEvolution, problem F = parameters.F CR = parameters.CR - X = zeros(N,D) + X = zeros(eltype(xBest), N,D) for i in 1:N x = get_position(population[i]) u = DE_mutation(population, F, strategy, 1) v = DE_crossover(x, u, CR) - evo_boundary_repairer!(v, xBest, problem.bounds) - X[i,:] = v + evo_boundary_repairer!(v, xBest, problem.search_space) + X[i,:] = _fix_type(v, problem.search_space) end X diff --git a/src/algorithms/DE/epsilonDE.jl b/src/algorithms/DE/epsilonDE.jl index b106a7ed..66846249 100644 --- a/src/algorithms/DE/epsilonDE.jl +++ b/src/algorithms/DE/epsilonDE.jl @@ -58,7 +58,7 @@ function initialize!( elite_sols = sortperm(status.population, lt = is_better) θ = round(Int, 0.2length(elite_sols)) - s = rand(status.population[elite_sols[1:θ]]) + s = rand(options.rng, status.population[elite_sols[1:θ]]) parameters.ε_0 = sum_violations(s) parameters.Tc = round(Int, 0.2*options.iterations) diff --git a/src/algorithms/ECA/CECA.jl b/src/algorithms/ECA/CECA.jl index 47c71725..9361287f 100644 --- a/src/algorithms/ECA/CECA.jl +++ b/src/algorithms/ECA/CECA.jl @@ -8,7 +8,7 @@ function update_state!( kargs... ) - I = randperm(parameters.N) + I = randperm(options.rng, parameters.N) N = parameters.N population = status.population @@ -22,7 +22,7 @@ function update_state!( feasible_solutions = findall( s->s.is_feasible, population ) weights = compute_weights(population) - X_next = zeros(parameters.N, size(problem.bounds, 2)) + X_next = zeros(parameters.N, getdim(problem)) # For each elements in Population for i = 1:parameters.N @@ -40,7 +40,7 @@ function update_state!( u_best = argmax(mass) # stepsize - η = parameters.η_max * rand() + η = parameters.η_max * rand(options.rng) # u: worst element in U u = U[u_worst].x @@ -49,10 +49,10 @@ function update_state!( # current-to-center/bin y = x .+ η .* (c .- u) - mask = rand(length(y)) .< 1.0 / length(y) + mask = rand(options.rng, length(y)) .< 1.0 / length(y) y[mask] = v[mask] - evo_boundary_repairer!(y, c, problem.bounds) + evo_boundary_repairer!(y, c, problem.search_space, options.rng) X_next[i,:] = y end diff --git a/src/algorithms/ECA/ECA.jl b/src/algorithms/ECA/ECA.jl index fcabe6cc..621cae31 100644 --- a/src/algorithms/ECA/ECA.jl +++ b/src/algorithms/ECA/ECA.jl @@ -69,8 +69,7 @@ function ECA(; ε::Float64 = 0.0, adaptive::Bool = false, resize_population::Bool = false, - information = Information(), - options = Options(), + kargs... ) @@ -89,22 +88,18 @@ function ECA(; adaptive, resize_population, ) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm(parameters; kargs...) end -function eca_solution(status, parameters, options, problem, I, i) +function eca_solution(status, parameters, options, problem, I, i, rng = default_rng_mh()) p = status.f_calls / options.f_calls_limit # generate U masses U = getU(status.population, parameters.K, I, i, parameters.N) # generate center of mass c, u_worst, u_best = center(U) # stepsize - η = parameters.η_max * rand() + η = parameters.η_max * rand(rng) # u: worst element in U u = U[u_worst] |> get_position # current-to-center/bin @@ -118,8 +113,8 @@ function eca_solution(status, parameters, options, problem, I, i) y = get_position(status.population[i]) .+ η .*(minimizer(status) .- c) end # binary crossover - y, M_current = crossover(U[u_best].x, y, parameters.p_cr) - evo_boundary_repairer!(y, c, problem.bounds) + y, M_current = crossover(U[u_best].x, y, parameters.p_cr, options.rng) + evo_boundary_repairer!(y, c, problem.search_space, rng) y, M_current end @@ -134,8 +129,8 @@ function update_state!( kargs... ) - I = randperm(parameters.N) - D = size(problem.bounds, 2) + I = randperm(options.rng, parameters.N) + D = getdim(problem) parameters.adaptive && (Mcr_fail = zeros(D)) @@ -145,7 +140,7 @@ function update_state!( # For each elements in Population for i in 1:parameters.N - y, M_current = eca_solution(status, parameters, options,problem,I,i) + y, M_current = eca_solution(status, parameters, options,problem,I,i, options.rng) if options.parallel_evaluation X_next[i,:] = y @@ -174,7 +169,7 @@ function update_state!( if parameters.adaptive parameters.p_cr = - adaptCrossover(parameters.p_cr, Mcr_fail / parameters.N) + adaptCrossover(parameters.p_cr, Mcr_fail / parameters.N, options.rng) end resize_population!(status, parameters, options) @@ -237,11 +232,18 @@ function initialize!( args...; kargs... ) - D = size(problem.bounds, 2) + + if !(problem.search_space isa BoxConstrainedSpace) + s = string(typeof(problem.search_space)) + error("ECA only works on box-contrained spaces but "*s*" was given.") + end + + + D = getdim(problem) if parameters.N <= parameters.K - parameters.N = parameters.K * D + parameters.N = min(parameters.K * D, 1000) end if options.f_calls_limit == 0 @@ -258,13 +260,12 @@ function initialize!( if parameters.adaptive - parameters.p_cr = rand(D) + parameters.p_cr = rand(options.rng, D) else parameters.p_cr = parameters.p_bin .* ones(D) end - # initialize!(problem, nothing, parameters, status, information, options) st = gen_initial_state(problem,parameters,information,options, status) st @@ -299,7 +300,7 @@ function reproduction(status, parameters::ECA, problem) parameters.K, parameters.η_max; i=i, - bounds = problem.bounds) + bounds = problem.search_space) end X diff --git a/src/algorithms/ECA/adaptive_parameters.jl b/src/algorithms/ECA/adaptive_parameters.jl index f1ba55a6..26d337e2 100644 --- a/src/algorithms/ECA/adaptive_parameters.jl +++ b/src/algorithms/ECA/adaptive_parameters.jl @@ -1,9 +1,9 @@ -function adaptCrossover(p_cr::Vector{Float64}, M::Vector{Float64}) +function adaptCrossover(p_cr::Vector{Float64}, M::Vector{Float64}, rng = default_rng_mh()) p_best = p_cr[argmin(M)] for i = 1:length(p_cr) if M[i] > 0.3 - pn = abs(p_best .+ 0.3 * randn()) + pn = abs(p_best .+ 0.3 * randn(rng)) if pn > 1.0 pn = 1.0 end diff --git a/src/algorithms/ECA/center_of_mass.jl b/src/algorithms/ECA/center_of_mass.jl index f724cb95..c7fb4d9d 100644 --- a/src/algorithms/ECA/center_of_mass.jl +++ b/src/algorithms/ECA/center_of_mass.jl @@ -1,8 +1,10 @@ ################################################ # Unconstrained ################################################ -function fitnessToMass(fitness::Vector{Float64}) - m = minimum(fitness) +function fitnessToMass(fitness::Vector{<:Real}) + m, M = extrema(fitness) + # avoid dividing by zero + m ≈ 0 && M ≈ 0 && return ones(eltype(fitness), length(fitness)) if m < 0 fitness = 2 * abs(m) .+ fitness @@ -18,22 +20,22 @@ end get the mass in a vector `m[i] = f[i] + 2*max_fitness*sum_vio[i] `, tolerance in equality constraints are given by `epsilon` """ -getMass(U::Array{xf_indiv,1}) = fitnessToMass(fvals(U)) +getMass(U::Array{xf_solution{Vector{T}},1}) where T <: Real = fitnessToMass(fvals(U)) -function getMass(U::Array{xfgh_indiv,1}) +function getMass(U::Array{xfgh_solution{Vector{T}},1}) where T <: Real fitness = fvals(U) - fitnessToMass(fitness + 2maximum(abs.(fitness))*sum_violations.(U)) + M = 2maximum(abs.(fitness)) + fitnessToMass(fitness + max(M, 100)*sum_violations.(U)) end +function getMass(U::Array{xFgh_solution{Vector{T}},1}) where T <: Real + fitness = sum(fvals(U), dims=2) |> vec + M = 2maximum(abs.(fitness)) + fitnessToMass(fitness + max(M, 100)*sum_violations.(U)) +end function center(U, mass) - d = length(U[1].x) - - c = zeros(Float64, d) - - for i = 1:length(mass) - c += mass[i] .* U[i].x - end + c = sum(m * get_position(u) for (m, u) in zip(mass, U)) return c / sum(mass) end @@ -73,7 +75,7 @@ end -function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int, feasible_solutions) +function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int, feasible_solutions, rng = default_rng_mh()) # at least half of the population is feasible to generate random centers if length(feasible_solutions) >= 0.5N || length(feasible_solutions) == 0 return getU_ids(K, I, i, N) @@ -88,20 +90,21 @@ function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int, feasible_solutions) U_ids = I[j.+1] end - push!(U_ids, rand(feasible_solutions)) + push!(U_ids, rand(rng, feasible_solutions)) return U_ids end function crossover( - x::Vector{Float64}, - y::Vector{Float64}, + x::Vector{<:Real}, + y::Vector{<:Real}, p_cr::Vector{Float64}, + rng = default_rng_mh() ) D = length(x) tmp2 = zeros(D) for j = 1:D - if rand() < p_cr[j] + if rand(rng) < p_cr[j] y[j] = x[j] tmp2[j] += 1 end @@ -118,10 +121,9 @@ Compute a solution using ECA variation operator, `K` is the number of solutions calculate the center of mass and `η_max` is the maximum stepsize. """ function ECA_operator( - # population::AbstractArray{xf_indiv}, K, η_max; - population, K, η_max; - i = rand(1:length(population)), - U = rand(population, K), + population, K, η_max, rng = default_rng_mh(); + i = rand(rng, 1:length(population)), + U = rand(rng, population, K), bounds = nothing ) @@ -131,7 +133,7 @@ function ECA_operator( c, u_worst, u_best = center(U) # stepsize - η = η_max * rand() + η = η_max * rand(rng) # u: worst element in U u = get_position(U[u_worst]) @@ -141,7 +143,7 @@ function ECA_operator( return y end - !isempty(bounds) && evo_boundary_repairer!(y, c, bounds) + evo_boundary_repairer!(y, c, bounds, rng) return y end diff --git a/src/algorithms/GA/GA.jl b/src/algorithms/GA/GA.jl index 9a721415..0ffb8339 100644 --- a/src/algorithms/GA/GA.jl +++ b/src/algorithms/GA/GA.jl @@ -8,6 +8,9 @@ mutable struct GA{T1, T2, T3, T4, T5} <: AbstractGA environmental_selection::T5 end +iscompatible(search_space::BitArraySpace, algorithm::GA) = true +iscompatible(search_space::PermutationSpace, algorithm::GA) = true + """ GA(; N = 100, @@ -41,7 +44,7 @@ f (generic function with 1 method) julia> dim = 10; -julia> optimize(f, repeat([false, true], 1, dim), GA()) +julia> optimize(f, BitArraySpace(dim), GA()) +=========== RESULT ==========+ iteration: 500 minimum: 0 @@ -62,7 +65,7 @@ julia> perm_size = 10; julia> ga = GA(;initializer = RandomPermutation(N=100), crossover=OrderCrossover(), mutation=SlightMutation()); -julia> optimize(f, zeros(Int,2,perm_size), ga) +julia> optimize(f, PermutationSpace(perm_size), ga) +=========== RESULT ==========+ iteration: 500 minimum: 0 @@ -135,8 +138,7 @@ function GA(; crossover = UniformCrossover(p = p_crossover), mutation = BitFlipMutation(p = p_mutation), environmental_selection = ElitistReplacement(), - options = Options(), - information = Information() + kargs... ) parameters = GA(initializer, @@ -146,13 +148,44 @@ function GA(; environmental_selection ) - Algorithm( - parameters, - information = information, - options = options - ) + Algorithm( parameters; kargs...) end +# this method is called in src/optimize.jl +function get_parameters( + f, + search_space::BoxConstrainedSpace, # real/integer encoding + ::Type{T} + ) where T <: GA + + N = 100 + GA(; + N = N, + initializer = RandomInBounds(;N), + selection = TournamentSelection(;N), + crossover = SBX(;bounds=search_space), + mutation = PolynomialMutation(;bounds=search_space), + environmental_selection = ElitistReplacement() + ) +end + + +function get_parameters( + f, + search_space::PermutationSpace, # permutation-based encoding + ::Type{T} + ) where T <: GA + + N = 100 + GA(; + N = N, + initializer = RandomPermutation(;N), + selection = TournamentSelection(;N), + crossover = OrderCrossover(), + mutation = SlightMutation(), + environmental_selection = ElitistReplacement() + ) +end function initialize!( status, diff --git a/src/algorithms/MCCGA/MCCGA.jl b/src/algorithms/MCCGA/MCCGA.jl index 0aafc47a..069ccaeb 100644 --- a/src/algorithms/MCCGA/MCCGA.jl +++ b/src/algorithms/MCCGA/MCCGA.jl @@ -89,17 +89,12 @@ function MCCGA(; maxsamples = 10_000, mutation = 1 / N, use_local_search = true, - information = Information(), - options = Options() + kargs... ) parameters = MCCGA(N, maxsamples, mutation, [], use_local_search) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm(parameters;kargs...) end function initialize!( @@ -120,10 +115,10 @@ function initialize!( options.f_calls_limit = options.iterations * parameters.N + 1 end - lower = problem.bounds[1,:] - upper = problem.bounds[2,:] + lower = problem.search_space.lb + upper = problem.search_space.ub - parameters.probvector = initialprobs(lower, upper, maxsamples = parameters.maxsamples) + parameters.probvector = initialprobs(lower, upper, options.rng, maxsamples = parameters.maxsamples) # sample vectors to create an initial State return gen_initial_state(problem,parameters,information,options,status) @@ -143,8 +138,8 @@ function update_state!( chsize = length(probvector) # parents - ch1 = sample(probvector) - ch2 = sample(probvector) + ch1 = sample(probvector, options.rng) + ch2 = sample(probvector, options.rng) # evaluate cost function sol1 = create_solution(floats(ch1), problem) @@ -203,7 +198,7 @@ function final_stage!( costfunction(x) = evaluate(x, problem) probvector = parameters.probvector - sampledvector = sample(probvector) + sampledvector = sample(probvector, options.rng) # initial solution for the local search initial_solution = floats(sampledvector) diff --git a/src/algorithms/MCCGA/utils.jl b/src/algorithms/MCCGA/utils.jl index 680f6b9e..e53654bd 100644 --- a/src/algorithms/MCCGA/utils.jl +++ b/src/algorithms/MCCGA/utils.jl @@ -4,24 +4,26 @@ ################################################################## function initialprobs( - lower::Array{T,1}, - upper::Array{T,1}; - maxsamples = 10000, -) where {T<:Number} + lower::Array{T,1}, + upper::Array{T,1}, + rng; + maxsamples = 10000, + ) where {T<:Number} + p = length(lower) @assert p == length(upper) bitlen = p * 32 mat = Array{Int8,2}(undef, maxsamples, bitlen) for tries = 1:maxsamples - randvalues = map((x, y) -> x + rand() * (y - x), lower, upper) + randvalues = map((x, y) -> x + rand(rng) * (y - x), lower, upper) mat[tries, :] .= bits(randvalues) end return colmeans(mat) end -function sample(probs::Array{T,1}) where {T<:Number} - return map(x -> if rand() < x +function sample(probs::Array{T,1}, rng) where {T<:Number} + return map(x -> if rand(rng) < x 1 else 0 diff --git a/src/algorithms/MOEAD_DE/MOEAD_DE.jl b/src/algorithms/MOEAD_DE/MOEAD_DE.jl index 88675509..663630a2 100644 --- a/src/algorithms/MOEAD_DE/MOEAD_DE.jl +++ b/src/algorithms/MOEAD_DE/MOEAD_DE.jl @@ -107,8 +107,7 @@ function MOEAD_DE(weights ; B = Array{Int}[], s1 = 0.01, s2 = 20.0, - information = Information(), - options = Options(), + kargs... ) @@ -131,15 +130,7 @@ function MOEAD_DE(weights ; initialize_closest_weight_vectors!(parameters) - alg = Algorithm( - parameters, - information = information, - options = options, - ) - - - - alg + Algorithm(parameters;kargs...) end @@ -166,7 +157,7 @@ function initialize!( status = gen_initial_state(problem,parameters,information,options,status) - D = size(problem.bounds, 2) + D = getdim(problem) parameters.nobjectives = length(status.population[1].f) parameters.p_m = parameters.p_m < 0.0 ? 1.0 / D : parameters.p_m @@ -193,7 +184,7 @@ function update_state!( F = parameters.F CR = parameters.CR - D = size(problem.bounds, 2) + D = getdim(problem) N = parameters.N @@ -211,7 +202,7 @@ function update_state!( # reproduction v = MOEAD_DE_reproduction(i, P_idx, population, parameters, problem) # repair - replace_with_random_in_bounds!(v, problem.bounds) + replace_with_random_in_bounds!(v, problem.search_space) h = create_solution(v, problem) # update z update_reference_point!(parameters.z, h) @@ -252,13 +243,13 @@ Perform Differential Evolution operators and polynomial mutation using three vec `a, b, c` and parameters `F, CR, p_m, η`, i.e., stepsize, crossover and mutation probability. """ -function MOEAD_DE_reproduction(a, b, c, F, CR, p_m, η, bounds) +function MOEAD_DE_reproduction(a, b, c, F, CR, p_m, η, bounds::BoxConstrainedSpace) D = length(a) # binomial crossover v = zeros(length(a)) - la = view(bounds, 1, :) - lb = view(bounds, 2, :) + la = bounds.lb + lb = bounds.ub # binomial crossover for j in 1:D @@ -309,7 +300,7 @@ function MOEAD_DE_reproduction(i, P_idx, population, parameters::MOEAD_DE, probl parameters.CR, parameters.p_m, parameters.η, - problem.bounds) + problem.search_space) end g_te_ap(gx, V, τ, s1, s2) = V < τ ? gx + s1*V^2 : gx + s1*τ^2 + s2*(V - τ) diff --git a/src/algorithms/NSGA2/NSGA2.jl b/src/algorithms/NSGA2/NSGA2.jl index 4bb9db33..3bea148d 100644 --- a/src/algorithms/NSGA2/NSGA2.jl +++ b/src/algorithms/NSGA2/NSGA2.jl @@ -70,16 +70,11 @@ function NSGA2(; p_cr = 0.9, η_m = 20, p_m = -1, - information = Information(), - options = Options(), + kargs... ) parameters = NSGA2(N, promote( Float64(η_cr), p_cr, η_m, p_m )...) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm( parameters; kargs...) end @@ -120,7 +115,7 @@ offspring. Return two vectors. """ function GA_reproduction(pa::AbstractVector{T}, pb::AbstractVector{T}, - bounds; + bounds::BoxConstrainedSpace; η_cr = 20, η_m = 15, p_cr = 0.9, @@ -156,7 +151,7 @@ Same that `GA_reproduction` but only returns one offspring. """ function GA_reproduction_half(pa::AbstractVector{T}, pb::AbstractVector{T}, - bounds; + bounds::BoxConstrainedSpace; η_cr = 20, η_m = 15, p_cr = 0.9, @@ -176,11 +171,22 @@ function GA_reproduction_half(pa::AbstractVector{T}, return c end +function GA_reproduction_half(pa::AbstractVector{T}, + pb::AbstractVector{T}, + bounds::AbstractMatrix; + kargs... + ) where T <: AbstractFloat + + b = BoxConstrainedSpace(lb = bounds[1,:], ub = bounds[2,:]) + GA_reproduction_half(pa, pb, b; kargs...) + +end + function reproduction(pa, pb, parameters::AbstractNSGA, problem) # crossover and mutation c1, c2 = GA_reproduction(get_position(pa), get_position(pb), - problem.bounds; + problem.search_space; η_cr = parameters.η_cr, p_cr = parameters.p_cr, η_m = parameters.η_m, @@ -203,7 +209,7 @@ function initialize!( args...; kargs... ) - D = size(problem.bounds, 2) + D = getdim(problem) if parameters.p_m < 0.0 parameters.p_m = 1.0 / D @@ -256,7 +262,7 @@ function reproduction(status, parameters::AbstractNSGA, problem) @assert !isempty(status.population) N_half = parameters.N - Q = zeros(2N_half, size(problem.bounds, 2)) + Q = zeros(2N_half, getdim(problem)) for i in 1:N_half pa = tournament_selection(status.population) @@ -264,7 +270,7 @@ function reproduction(status, parameters::AbstractNSGA, problem) c1, c2 = GA_reproduction(get_position(pa), get_position(pb), - problem.bounds; + problem.search_space; η_cr = parameters.η_cr, p_cr = parameters.p_cr, η_m = parameters.η_m, diff --git a/src/algorithms/NSGA3/NSGA3.jl b/src/algorithms/NSGA3/NSGA3.jl index 14ba9dcb..cc8195fc 100644 --- a/src/algorithms/NSGA3/NSGA3.jl +++ b/src/algorithms/NSGA3/NSGA3.jl @@ -68,17 +68,12 @@ function NSGA3(; p_m = -1, partitions=12, reference_points=Vector{Float64}[], - information = Information(), - options = Options(), + kargs... ) parameters = NSGA3(N, promote( Float64(η_cr), p_cr, η_m, p_m )..., partitions, reference_points) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm(parameters; kargs...) end @@ -272,8 +267,9 @@ function initialize!( options::Options, args...; kargs... -) - D = size(problem.bounds, 2) + ) + + D = getdim(problem) if parameters.p_m < 0.0 parameters.p_m = 1.0 / D @@ -324,7 +320,7 @@ function reproduction(status, parameters::NSGA3, problem) @assert !isempty(status.population) I = randperm(parameters.N) - Q = zeros(parameters.N, size(problem.bounds, 2)) + Q = zeros(parameters.N, getdim(problem)) for i = 1:parameters.N ÷ 2 pa = status.population[I[2i-1]] @@ -332,7 +328,7 @@ function reproduction(status, parameters::NSGA3, problem) c1, c2 = GA_reproduction(get_position(pa), get_position(pb), - problem.bounds; + problem.search_space; η_cr = parameters.η_cr, p_cr = parameters.p_cr, η_m = parameters.η_m, diff --git a/src/algorithms/PSO/PSO.jl b/src/algorithms/PSO/PSO.jl index eb115247..b44e61d6 100644 --- a/src/algorithms/PSO/PSO.jl +++ b/src/algorithms/PSO/PSO.jl @@ -50,23 +50,18 @@ julia> optimize(f, [-1 -1 -1; 1 1 1.0], PSO(N = 100, C1=1.5, C2=1.5, ω = 0.7)) """ function PSO(; - N::Int = 0, - C1 = 2.0, - C2 = 2.0, - ω = 0.8, - v = Float64[], - flock = xf_indiv[], - information = Information(), - options = Options(), -) - -parameters = PSO(N, promote(Float64(C1), C2, ω)..., v, flock) - -Algorithm( - parameters, - information = information, - options = options, -) + N::Int = 0, + C1 = 2.0, + C2 = 2.0, + ω = 0.8, + v = Float64[], + flock = xf_indiv[], + kargs... + ) + + parameters = PSO(N, promote(Float64(C1), C2, ω)..., v, flock) + + Algorithm( parameters; kargs...) end function update_state!( @@ -80,15 +75,15 @@ function update_state!( ) xGBest = get_position(status.best_sol) - X_new = zeros(parameters.N, size(problem.bounds, 2)) + X_new = zeros(parameters.N, getdim(problem)) # For each elements in population for i in 1:parameters.N x = get_position(parameters.flock[i]) xPBest = get_position(status.population[i]) - parameters.v[i, :] = velocity(x, parameters.v[i, :], xPBest, xGBest, parameters) + parameters.v[i, :] = velocity(x, parameters.v[i, :], xPBest, xGBest, parameters, options.rng) x += parameters.v[i, :] - reset_to_violated_bounds!(x, problem.bounds) + reset_to_violated_bounds!(x, problem.search_space) X_new[i,:] = x end @@ -119,7 +114,7 @@ function initialize!( args...; kargs... ) - D = size(problem.bounds, 2) + D = getdim(problem) if parameters.N == 0 diff --git a/src/algorithms/PSO/velocity.jl b/src/algorithms/PSO/velocity.jl index 71fb065b..013d7d25 100644 --- a/src/algorithms/PSO/velocity.jl +++ b/src/algorithms/PSO/velocity.jl @@ -1,10 +1,10 @@ -function PSO_velocity(x, v, pbest, gbest, C1, C2, ω) - r1 = C1 * rand() - r2 = C2 * rand() +function PSO_velocity(x, v, pbest, gbest, C1, C2, ω, rng) + r1 = C1 * rand(rng) + r2 = C2 * rand(rng) ω * v + r1 * (pbest - x) + r2 * (gbest - x) end -function velocity(x, v, pbest, gbest, parameters) - PSO_velocity(x, v, pbest, gbest, parameters.C1, parameters.C2, parameters.ω) +function velocity(x, v, pbest, gbest, parameters, rng) + PSO_velocity(x, v, pbest, gbest, parameters.C1, parameters.C2, parameters.ω, rng) end diff --git a/src/algorithms/SA/SA.jl b/src/algorithms/SA/SA.jl index 9c57de3e..2ad36d40 100644 --- a/src/algorithms/SA/SA.jl +++ b/src/algorithms/SA/SA.jl @@ -65,17 +65,12 @@ function SA(; x_initial::Vector = zeros(0), N::Int = 500, tol_fun::Real= 1e-4, - information = Information(), - options = Options() + kargs... ) parameters = SA(N, x_initial, tol_fun, zeros(0), Inf) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm( parameters; kargs...) end @@ -91,11 +86,12 @@ function initialize!( ) - l = view( problem.bounds, 1,:) - u = view( problem.bounds, 2,:) + l = problem.search_space.lb + u = problem.search_space.ub + rng = options.rng if isempty(parameters.x_initial) - parameters.x_initial = l .+ (u .- l) .* rand(length(u)) + parameters.x_initial = l .+ (u .- l) .* rand(rng, length(u)) end if options.f_calls_limit == 0 @@ -134,26 +130,27 @@ function update_state!( max_evals = options.f_calls_limit N = parameters.N TolFun = parameters.tol_fun + rng = options.rng # T is the inverse of temperature. T = nevals / max_evals μ = 10.0 ^( 100T ) - l = view( problem.bounds, 1,:) - u = view( problem.bounds, 2,:) + l = problem.search_space.lb + u = problem.search_space.ub # For each temperature we take 500 test points to simulate reach termal # equilibrium. for i = 1:N # We generate new test point using newSol function - dx = newSol(2rand(length(parameters.x)) .- 1.0 , μ) .* (u-l) + dx = newSol(2rand(rng, length(parameters.x)) .- 1.0 , μ) .* (u-l) # the test point and fx1=f(x1) x1 = parameters.x + dx # Next step is to keep solution within bounds #x1 = (x1 .< l).*l+(l .<= x1).*(x1 .<= u).*x1+(u .< x1).*u - reset_to_violated_bounds!(x1, problem.bounds) + reset_to_violated_bounds!(x1, problem.search_space) fx1 = evaluate(x1, problem) status.f_calls += 1 @@ -164,7 +161,7 @@ function update_state!( # If the function variation,df, is <0 we take test point as current # point. And if df>0 we use Metropolis condition to accept or # reject the test point as current point. - if (df < 0 || rand() < exp(-T*df/(abs(parameters.fx)) / TolFun)) + if (df < 0 || rand(rng) < exp(-T*df/(abs(parameters.fx)) / TolFun)) parameters.x = x1 parameters.fx= fx1 end diff --git a/src/algorithms/SMS_EMOA/SMS_EMOA.jl b/src/algorithms/SMS_EMOA/SMS_EMOA.jl index 5f8eafcd..174fa006 100644 --- a/src/algorithms/SMS_EMOA/SMS_EMOA.jl +++ b/src/algorithms/SMS_EMOA/SMS_EMOA.jl @@ -72,16 +72,11 @@ function SMS_EMOA(; η_m = 20, p_m = -1, n_samples = 10_000, - information = Information(), - options = Options(), + kargs... ) parameters = SMS_EMOA(N, promote( Float64(η_cr), p_cr, η_m, p_m)..., n_samples) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm(parameters; kargs...) end @@ -102,7 +97,7 @@ function update_state!( J = randperm(parameters.N) if options.parallel_evaluation - Q = zeros(parameters.N, size(problem.bounds, 2)) + Q = zeros(parameters.N, getdim(problem)) end for i = 1:parameters.N @@ -112,15 +107,15 @@ function update_state!( # crossover _, c = SBX_crossover(get_position(pa), get_position(pb), - problem.bounds, + problem.search_space, parameters.η_cr, parameters.p_cr) # mutation - polynomial_mutation!(c, problem.bounds, parameters.η_m, parameters.p_m) + polynomial_mutation!(c, problem.search_space, parameters.η_m, parameters.p_m) # repair solutions if necesary - reset_to_violated_bounds!(c, problem.bounds) + reset_to_violated_bounds!(c, problem.search_space) if options.parallel_evaluation Q[i,:] = c @@ -165,9 +160,9 @@ function initialize!( information::Information, options::Options, args...; - kargs... -) - D = size(problem.bounds, 2) + kargs...) + + D = getdim(problem) if parameters.p_m < 0.0 parameters.p_m = 1.0 / D diff --git a/src/algorithms/SPEA2/SPEA2.jl b/src/algorithms/SPEA2/SPEA2.jl index 80706414..169ebe5c 100644 --- a/src/algorithms/SPEA2/SPEA2.jl +++ b/src/algorithms/SPEA2/SPEA2.jl @@ -67,17 +67,11 @@ function SPEA2(; p_cr = 0.9, η_m = 20, p_m = -1, - information = Information(), - options = Options(), + kargs... ) parameters = SPEA2(N, promote( Float64(η_cr), p_cr, η_m, p_m )...,[]) - Algorithm( - parameters, - information = information, - options = options, - ) - + Algorithm(parameters; kargs...) end diff --git a/src/algorithms/WOA/WOA.jl b/src/algorithms/WOA/WOA.jl index 3b826ef1..2bfe96f6 100644 --- a/src/algorithms/WOA/WOA.jl +++ b/src/algorithms/WOA/WOA.jl @@ -51,15 +51,11 @@ julia> optimize(f, [-1 -1 -1; 1 1 1.0], WOA(N = 100)) ``` """ -function WOA(;N = 30, information = Information(), options = Options()) +function WOA(;N = 30, kargs...) parameters = WOA(N) - Algorithm( - parameters, - information = information, - options = options, - ) + Algorithm( parameters; kargs...) end @@ -73,7 +69,8 @@ function initialize!( kargs... ) - lb, ub= problem.bounds[1,:], problem.bounds[2,:] + lb = problem.search_space.lb + ub = problem.search_space.ub #Initialize the positions of search agents max_it = 500 @@ -110,8 +107,9 @@ function update_state!( Max_iter = options.iterations N = parameters.N - D = size(problem.bounds, 2) + D = getdim(problem) t = status.iteration + rng = options.rng a=2-t*((2)/Max_iter) # a decreases linearly fron 2 to 0 in Eq. (2.3) @@ -122,23 +120,23 @@ function update_state!( # Update the Position of search agents for i=1:N - r1 = rand() # r1 is a random number in [0,1] - r2 = rand() # r2 is a random number in [0,1] + r1 = rand(rng) # r1 is a random number in [0,1] + r2 = rand(rng) # r2 is a random number in [0,1] A = 2a*r1 - a # Eq. (2.3) in the paper C = 2r2 # Eq. (2.4) in the paper b = 1 # parameters in Eq. (2.5) - l = (a2-1)*rand()+1 # parameters in Eq. (2.5) + l = (a2-1)*rand(rng)+1 # parameters in Eq. (2.5) - p = rand() # p in Eq. (2.6) + p = rand(rng) # p in Eq. (2.6) - x = copy(status.population[i].x) + x = zeros(D) for j = 1:D if p < 0.5 if abs(A) >= 1 - rand_leader_index = floor(Int, N* rand()+1) + rand_leader_index = floor(Int, N* rand(rng)+1) X_rand = status.population[rand_leader_index].x D_X_rand= abs(C*X_rand[j] - x[j] ) # Eq. (2.7) x[j] = X_rand[j]-A*D_X_rand # Eq. (2.8) @@ -158,7 +156,7 @@ function update_state!( end # for j - reset_to_violated_bounds!(x, problem.bounds) + reset_to_violated_bounds!(x, problem.search_space) X_new[i,:] = x end # for i diff --git a/src/algorithms/algorithm.jl b/src/algorithms/algorithm.jl index e69de29b..ec60392f 100644 --- a/src/algorithms/algorithm.jl +++ b/src/algorithms/algorithm.jl @@ -0,0 +1,86 @@ +mutable struct Algorithm{T} <: AbstractAlgorithm + parameters::T + status::State + information::Information + options::Options + termination::Termination +end + +function Algorithm( + parameters; + initial_state::State = State(nothing, []), + information::Information = Information(), + options::Options = Options(), + termination = Termination(), + kargs... + ) + if !isempty(kargs) + _pp = collect(kargs) + p = string.(first.(_pp)) .* " = " .* string.(last.(_pp)) + pp = join(p, ", ") + @warn "Parameters: `$pp` never used." + end + + if !(termination isa Termination) + termination = Termination(checkany = [termination]) + end + + + Algorithm(parameters, initial_state, information, options, termination) + +end + +function _print_title(io, title) + decor = join(repeat('=', length(title))) + printstyled(io, title, "\n", color=:blue, bold=true) + printstyled(io, decor, "\n", color=:gray) +end + +function Base.show(io::IO, alg::Algorithm) + _print_title(io, "Algorithm Parameters") + print(io, " ") + Base.show(IOContext(io, :compact => true), alg.parameters) + println(io, "\n") + + show(io, alg.status) +end + +function Base.show(io::IO, parameters::AbstractParameters) + s = typeof(parameters) + + fns = fieldnames(s) + all_vals = [getfield(parameters, f) for f in fns] + # remove those parameters containing empty arrays before showing + mask = findall(v -> !(v isa Array), all_vals) + vals = [sprint(show, v) for v in all_vals[mask]] + str = string(s) * "(" * join(string.(fns[mask]) .* "=" .* vals, ", ") * ")" + + print(io, str) +end + +termination_status_message(alg::Algorithm) = termination_status_message(alg.status) +should_stop(algorithm::AbstractAlgorithm) = algorithm.status.stop +function get_result(algorithm::AbstractAlgorithm) + if isnothing(algorithm.status.best_sol) + error("First optimize a function: `optimize!(f, bounds, method)`") + end + + algorithm.status +end + + + +iscompatible(search_space, algorithm) = false +iscompatible(search_space::BoxConstrainedSpace, algorithm::AbstractParameters) = true +iscompatible(search_space::AbstractMatrix, algorithm::AbstractParameters) = true + +function check_compat(search_space, algorithm::AbstractParameters) + if iscompatible(search_space, algorithm) + return + end + + a = typeof(algorithm) + s = typeof(search_space) + error("Implemented metaheuristic `$a` is not compatible with search space `$s`. Try using bounds.") +end + diff --git a/src/algorithms/stop_criteria.jl b/src/algorithms/stop_criteria.jl index ce13c3ab..93128a1a 100644 --- a/src/algorithms/stop_criteria.jl +++ b/src/algorithms/stop_criteria.jl @@ -1,3 +1,4 @@ +#= function stop_criteria!( status, parameters::Union{ECA, DE, PSO, CGSA}, @@ -11,3 +12,4 @@ function stop_criteria!( status.stop = diff_check(status, information, options) end +=# diff --git a/src/common/gen_initial_state.jl b/src/common/gen_initial_state.jl index 7404b213..16de5236 100644 --- a/src/common/gen_initial_state.jl +++ b/src/common/gen_initial_state.jl @@ -13,7 +13,7 @@ function gen_initial_state(problem,parameters,information,options, status) options.debug && @warn("Population size in provided State differs from that in parameters") end - size(problem.bounds,2) != length(get_position(status.best_sol)) && + getdim(problem) != length(get_position(status.best_sol)) && error("Invalid population (dimension does not match with bounds)") _complete_population!(status,problem,parameters,information,options) @@ -29,7 +29,7 @@ end function gen_initial_state(problem,parameters,information,options) # population array - population = generate_population(parameters.N, problem,ε=options.h_tol) + population = generate_population(parameters.N, problem, options.rng,ε=options.h_tol) # best solution best_solution = get_best(population) @@ -37,14 +37,6 @@ function gen_initial_state(problem,parameters,information,options) return State(best_solution, population; f_calls = length(population), iteration=1) end -function Base.show(io::IO, parameters::AbstractParameters) - s = typeof(parameters) - - vals = string.(map(f -> getfield(parameters, f), fieldnames(s))) - str = string(s) * "(" * join(string.(fieldnames(s)) .* "=" .* vals, ", ") * ")" - - print(io, str) -end function _complete_population!(status,problem,parameters,information,options) if parameters.N < length(status.population) diff --git a/src/common/repair.jl b/src/common/repair.jl index 2e1c31e3..f2c1ee5b 100644 --- a/src/common/repair.jl +++ b/src/common/repair.jl @@ -1,7 +1,9 @@ -function reflected_back_to_bound!(x, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] +function reflected_back_to_bound!(x, bounds::BoxConstrainedSpace) + !bounds.rigid && (return x) + + for i in 1:getdim(bounds) + l = bounds.lb[i] + u = bounds.ub[i] while !( l <= x[i] <= u ) x[i] = x[i] < l ? 2l - x[i] : 2u - x[i] end @@ -10,22 +12,31 @@ function reflected_back_to_bound!(x, bounds) x end -function replace_with_random_in_bounds!(x, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] - if !( l <= x[i] <= u ) - x[i] = l + rand() * (u - l) - end +reflected_back_to_bound!(x, bounds::AbstractMatrix) = reflected_back_to_bound!(x, _mat_to_bounds(bounds)) + +function replace_with_random_in_bounds!(x, bounds::BoxConstrainedSpace, rng = default_rng_mh()) + !bounds.rigid && (return x) + lb = bounds.lb + ub = bounds.ub + Δ = bounds.Δ + for i in 1:getdim(bounds) + if !( lb[i] <= x[i] <= ub[i] ) + x[i] = lb[i] + rand(rng) * Δ[i] + end end x end -function wrap_to_bounds!(x, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] +replace_with_random_in_bounds!(x, bounds::AbstractMatrix) = replace_with_random_in_bounds!(x, _mat_to_bounds(bounds)) + +function wrap_to_bounds!(x, bounds::BoxConstrainedSpace) + !bounds.rigid && (return x) + lb = bounds.lb + ub = bounds.ub + for i in 1:getdim(bounds) + l = lb[i] + u = ub[i] if !( l <= x[i] <= u ) ρ = u - l x[i] = x[i] < l ? u - (l - x[i]) % ρ : l + (x[i] - u) % ρ @@ -34,11 +45,15 @@ function wrap_to_bounds!(x, bounds) end x end +wrap_to_bounds!(x, bounds::AbstractMatrix) = wrap_to_bounds!(x, _mat_to_bounds(bounds)) -function reset_to_violated_bounds!(x, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] +function reset_to_violated_bounds!(x, bounds::BoxConstrainedSpace) + !bounds.rigid && (return x) + lb = bounds.lb + ub = bounds.ub + for i in 1:getdim(bounds) + l = lb[i] + u = ub[i] if l > x[i] x[i] = l elseif x[i] > u @@ -48,28 +63,33 @@ function reset_to_violated_bounds!(x, bounds) end x end +reset_to_violated_bounds!(x, bounds::AbstractMatrix) = reset_to_violated_bounds!(x, _mat_to_bounds(bounds)) -function evo_boundary_repairer!(x, x_best, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] +function evo_boundary_repairer!(x, x_best, bounds::BoxConstrainedSpace, rng=default_rng_mh()) + !bounds.rigid && (return x) + lb = bounds.lb + ub = bounds.ub + for i in 1:getdim(bounds) + l = lb[i] + u = ub[i] if l > x[i] - α = rand() + α = rand(rng) x[i] = α*l + (1.0 - α) * x_best[i] elseif x[i] > u - β = rand() + β = rand(rng) x[i] = β*u + (1.0 - β) * x_best[i] end end x end +evo_boundary_repairer!(x, x_best, bounds::AbstractMatrix, rng=default_rng_mh()) = evo_boundary_repairer!(x, x_best, _mat_to_bounds(bounds), rng) -function is_in_bounds(x, bounds) - for i in 1:size(bounds,2) - l = bounds[1, i] - u = bounds[2, i] +function is_in_bounds(x, bounds::BoxConstrainedSpace) + for i in 1:getdim(bounds) + l = bounds.lb[i] + u = bounds.ub[i] if !( l <= x[i] <= u ) return false end @@ -78,3 +98,4 @@ function is_in_bounds(x, bounds) return true end +is_in_bounds(x, bounds::AbstractMatrix) = is_in_bounds(x, _mat_to_bounds(bounds)) diff --git a/src/common/show_plots.jl b/src/common/show_plots.jl index de15d5aa..b95e0948 100644 --- a/src/common/show_plots.jl +++ b/src/common/show_plots.jl @@ -86,3 +86,10 @@ end Base.show(io::IO, ::MIME"text/html", pop::Vector{Bee{xf_indiv}}) = show(io, "text/html", [p.sol for p in pop]) + + +function show_pf(io, pf::Vector{T}) where T <: xFgh_solution + show(io, "text/plain", pf) + print(io, "\n") +end + diff --git a/src/common/utils.jl b/src/common/utils.jl index 523b7412..a5518ce3 100644 --- a/src/common/utils.jl +++ b/src/common/utils.jl @@ -10,3 +10,21 @@ function pairwise_distances(A::Matrix, dist = Euclidean(); diag_val=Inf) return D end +function _mat_to_bounds(bounds::AbstractMatrix) + # validate bounds in cols + if size(bounds, 1) > 2 && size(bounds, 2) == 2 + bounds = bounds' + elseif size(bounds, 1) != 2 + error("Provide valid bounds. Suggestion: set `bounds = BoxConstrainedSpace(lb, ub)`.") + end + + lb = bounds[1,:] + ub = bounds[2,:] + BoxConstrainedSpace(lb, ub) +end +_mat_to_bounds(space::AbstractSearchSpace) = space +function _mat_to_bounds(b::Tuple{AbstractVector, AbstractVector}) + @assert length(b) == 2 "Provide valid lower and upper bounds." + BoxConstrainedSpace(first(b), last(b)) +end + diff --git a/src/core/abstracts.jl b/src/core/abstracts.jl index 803a4938..19e40cf2 100644 --- a/src/core/abstracts.jl +++ b/src/core/abstracts.jl @@ -2,4 +2,23 @@ abstract type AbstractSolution end abstract type AbstractAlgorithm end abstract type AbstractParameters end abstract type AbstractProblem end +abstract type AbstractTermination end + +""" + TerminationStatusCode +An Enum of possible => values for [`State`](@ref). +Possible values: + +- `ITERATION_LIMIT` +- `TIME_LIMIT` +- `EVALUATIONS_LIMIT` +- `ACCURACY_LIMIT` +- `OBJECTIVE_VARIANCE_LIMIT` +- `OBJECTIVE_DIFFERENCE_LIMIT` +- `OTHER_LIMIT` +- `UNKNOWN_STOP_REASON` + +See also [`termination_status_message`](@ref). +""" +const TerminationStatusCode = AbstractTermination diff --git a/src/core/options.jl b/src/core/options.jl index 6a64cc88..ca25b7c0 100644 --- a/src/core/options.jl +++ b/src/core/options.jl @@ -1,3 +1,5 @@ +default_rng_mh(seed=1223334444) = Random.default_rng(Int(seed)) + ##################################################### # # STRUCTURES FOR THE OPTIONS @@ -8,35 +10,35 @@ mutable struct Options f_tol::Float64 g_tol::Float64 h_tol::Float64 + f_tol_rel::Float64 f_calls_limit::Float64 - g_calls_limit::Float64 - h_calls_limit::Float64 time_limit::Float64 iterations::Int store_convergence::Bool - show_results::Bool debug::Bool - search_type::Symbol seed::UInt parallel_evaluation::Bool + verbose::Bool + rng end """ Options(; - x_tol::Real = 0.0, - f_tol::Real = 0.0, + x_tol::Real = 1e-8, + f_tol::Real = 1e-12, + f_tol_rel::Real = eps(), + f_tol_abs::Real = 0.0, g_tol::Real = 0.0, h_tol::Real = 0.0, f_calls_limit::Real = 0, - g_calls_limit::Real = 0, - h_calls_limit::Real = 0, time_limit::Real = Inf, iterations::Int = 1, store_convergence::Bool = false, debug::Bool = false, seed = rand(UInt), parallel_evaluation = false, + verbose = false, ) @@ -47,18 +49,21 @@ Main properties: - `x_tol` tolerance to the true minimizer if specified in `Information`. - `f_tol` tolerance to the true minimum if specified in `Information`. +- `f_tol_rel` relative tolerance. - `f_calls_limit` is the maximum number of function evaluations limit. - `time_limit` is the maximum time that `optimize` can spend in seconds. -- `iterations` is the maximum number iterationn permited. +- `iterations` is the maximum number of allowed iterations. - `store_convergence` if `true`, then push the current `State` in `State.convergence` at each generation/iteration - `debug` if `true`, then `optimize` function reports the current `State` (and interest information) for each iterations. - `seed` non-negative integer for the random generator seed. +- `parallel_evaluation` enables batch evaluations. +- `verbose` show simplified results each iteration. +- `rng` user-defined Random Number Generator. # Example -```jldoctest -julia> options = Options(f_calls_limit = 1000, debug=true, seed=1) -Options(0.0, 0.0, 0.0, 0.0, 1000.0, 0.0, 0.0, 0, false, true, true, :minimize, 0x0000000000000001) +```julia-repl +julia> options = Options(f_calls_limit = 1000, debug=false, seed=1); julia> f(x) = sum(x.^2) f (generic function with 1 method) @@ -69,35 +74,17 @@ julia> bounds = [ -10.0 -10 -10; # lower bounds -10.0 -10.0 -10.0 10.0 10.0 10.0 -julia> state = optimize(f, bounds, ECA(options=options)) -[ Info: Initializing population... -[ Info: Starting main loop... -+=========== RESULT ==========+ -| Iter.: 1 -| f(x) = 6.97287 -| solution.x = [-2.3628796262231875, -0.6781207370770752, -0.9642728360479853] -| f calls: 42 -| Total time: 0.0004 s -+============================+ - -... - -[ Info: Stopped since call_limit was met. -+=========== RESULT ==========+ -| Iter.: 47 -| f(x) = 1.56768e-08 -| solution.x = [-2.2626761322304715e-5, -9.838697194048792e-5, 7.405966506272336e-5] -| f calls: 1000 -| Total time: 0.0313 s -+============================+ +julia> state = optimize(f, bounds, ECA(options=options)); ``` """ function Options(; - x_tol = 0.0, - f_tol = 0.0, + x_tol = 1e-8, + f_tol = 1e-12, g_tol = 0.0, h_tol = 0.0, + f_tol_rel = eps(), + f_tol_abs = 0.0, f_calls_limit = 0.0, g_calls_limit = 0.0, h_calls_limit = 0.0, @@ -108,20 +95,37 @@ function Options(; debug = false, search_type = :minimize, parallel_evaluation = false, - seed = rand(UInt) + verbose = false, + seed = rand(UInt32), + rng = default_rng_mh(seed), ) Options( - promote(x_tol, f_tol, g_tol, h_tol)..., - promote(f_calls_limit, g_calls_limit, h_calls_limit, time_limit)..., + promote(x_tol, f_tol, g_tol, h_tol, f_tol_rel)..., + promote(f_calls_limit, time_limit)..., promote(iterations)..., # Results options - promote(store_convergence, show_results, debug)..., - Symbol(search_type), + promote(store_convergence, debug)..., UInt(seed), - Bool(parallel_evaluation) + Bool(parallel_evaluation), + Bool(verbose), + rng ) end + +function Base.show(io::IO, options::Options) + _print_title(io, "Options") + + txt = String[] + ln = Int[] + for field in fieldnames(Options) + v = getfield(options, field) + fl = string(field) .* ":" + push!(txt, @sprintf(" %-16s %s\n", fl, string(v))) + push!(ln, length(fl)) + end + println(io, join(txt[sortperm(ln)])) +end diff --git a/src/core/state.jl b/src/core/state.jl index 6b8c72f8..5fcc5077 100644 --- a/src/core/state.jl +++ b/src/core/state.jl @@ -104,13 +104,13 @@ end minimizer(state) Returns the approximation to the minimizer (argmin f(x)) stored in `state`. """ -minimizer(s::State) = s.best_sol.x +minimizer(s::State) = get_position(s.best_sol) """ minimum(state::Metaheuristics.State) Returns the approximation to the minimum (min f(x)) stored in `state`. """ -minimum(s::State) = s.best_sol.f +minimum(s::State) = fval(s.best_sol) """ positions(state) @@ -161,7 +161,7 @@ julia> n_fes, fxs = convergence(state); ``` """ -convergence(s::State) = begin +function convergence(s::State) sc = s.convergence isempty(sc) ? (zeros(Int, 0), zeros(0)) : (map(nfes, sc), map(minimum, sc)) end diff --git a/src/core/stop_status_codes.jl b/src/core/stop_status_codes.jl index b04a3b02..17b4d396 100644 --- a/src/core/stop_status_codes.jl +++ b/src/core/stop_status_codes.jl @@ -14,16 +14,26 @@ Possible values: See also [`termination_status_message`](@ref). """ -@enum(TerminationStatusCode, - ITERATION_LIMIT, - TIME_LIMIT, - EVALUATIONS_LIMIT, - ACCURACY_LIMIT, - OBJECTIVE_VARIANCE_LIMIT, - OBJECTIVE_DIFFERENCE_LIMIT, - OTHER_LIMIT, - UNKNOWN_STOP_REASON, - ) +const TerminationStatusCode = AbstractTermination + +struct IterationLimit <: TerminationStatusCode end +struct TimeLimit <: TerminationStatusCode end +struct EvaluationsLimit <: TerminationStatusCode end +struct AccuracyLimit <: TerminationStatusCode end +struct OtherLimit <: TerminationStatusCode end +struct UnknownStopReason <: TerminationStatusCode end +struct ObjectiveVarianceLimit <: TerminationStatusCode end +struct ObjectiveDifferenceLimit <: TerminationStatusCode end + + +const ITERATION_LIMIT = IterationLimit() +const TIME_LIMIT = TimeLimit() +const EVALUATIONS_LIMIT = EvaluationsLimit() +const ACCURACY_LIMIT = AccuracyLimit() +const OTHER_LIMIT = OtherLimit() +const UNKNOWN_STOP_REASON = UnknownStopReason() +const OBJECTIVE_VARIANCE_LIMIT = ObjectiveVarianceLimit() +const OBJECTIVE_DIFFERENCE_LIMIT= ObjectiveDifferenceLimit() """ termination_status_message(status) diff --git a/src/core/structures.jl b/src/core/structures.jl index f22c4022..deb2b7c9 100644 --- a/src/core/structures.jl +++ b/src/core/structures.jl @@ -1,52 +1,49 @@ -abstract type AbstractAlgorithm end +# Base.@kwdef struct Termination +# criteria_all::Vector = Any[] +# criteria_any::Vector = Any[] +# end -mutable struct Algorithm{T} <: AbstractAlgorithm - parameters::T - status::State - information::Information - options::Options -end - -function Algorithm( - parameters; - initial_state::State = State(nothing, []), - information::Information = Information(), - options::Options = Options(), -) - - Algorithm(parameters, initial_state, information, options) +mutable struct Problem{S} <: AbstractProblem + f::Function + search_space::S + f_calls::Int + parallel_evaluation::Bool end -function Base.show(io::IO, alg::Algorithm) - Base.show(io, alg.parameters) +function Base.getproperty(obj::Problem, sym::Symbol) + if sym === :bounds + @warn "`bounds` property is deprecated. Use `search_space` instead" + se = obj.search_space + return Array([se.lb se.ub]') + else # fallback to getfield + return getfield(obj, sym) + end end -termination_status_message(alg::Algorithm) = termination_status_message(alg.status) -should_stop(algorithm::AbstractAlgorithm) = algorithm.status.stop -function get_result(algorithm::AbstractAlgorithm) - if isnothing(algorithm.status.best_sol) - error("First optimize a function: `optimize!(f, bounds, method)`") - end - - algorithm.status +function getdim(problem::Problem) + return SearchSpaces.getdim(problem.search_space) end -mutable struct Problem{T} <: AbstractProblem - f::Function - bounds::Array{T,2} - f_calls::Int - parallel_evaluation::Bool + +function Problem(f::Function, search_space::AbstractSearchSpace; parallel_evaluation=false) + Problem(f, search_space, 0, parallel_evaluation) end -function Problem(f::Function, bounds::Array; parallel_evaluation=false) +function Problem(f::Function, bounds::AbstractMatrix{T}; kargs...) where T <: Number + # old problem definition if size(bounds,1) > 2 && size(bounds,2) == 2 bounds = Array(bounds') end - Problem(f, bounds, 0, parallel_evaluation) + Problem(f, BoxConstrainedSpace(lb = bounds[1,:], ub = bounds[2,:]); kargs...) +end + + +function Problem(f::Function, bounds::Array{Bool,2}; kargs...) + Problem(f, BitArraySpace(size(bounds,2)); kargs...) end function evaluate(x, problem::Problem) @@ -54,4 +51,3 @@ function evaluate(x, problem::Problem) return problem.f(x) end - diff --git a/src/externals.jl b/src/externals.jl index 897d6afa..f481863a 100644 --- a/src/externals.jl +++ b/src/externals.jl @@ -1,5 +1,6 @@ +import Random import Random: randperm, shuffle!, shuffle, seed! -import Printf.@printf +import Printf: @printf, @sprintf import LinearAlgebra import LinearAlgebra: norm, Diagonal, dot import Statistics: var, mean, std @@ -8,6 +9,29 @@ import Base.minimum using JMcDM using Requires +using Reexport + +@reexport using SearchSpaces +import SearchSpaces: AbstractSearchSpace,getdim +import SearchSpaces: AtomicSearchSpace, AbstractSampler, Sampler + +""" + boxconstraints(lb, ub, [rigid]) + BoxConstrainedSpace(ub, lb, [rigid]) + +Define a box-constrained search space (alias for BoxConstrainedSpace). + +See [`SearchSpaces.BoxConstrainedSpace`](@ref) for more details. + +# Example + +```julia +f(x) = (x[1] - 100)^2 + sum(abs.(x[2:end])) +bounds = boxconstraints(lb = zeros(5), ub = ones(5), rigid = false) +optimize(f, bounds, ECA) +``` +""" +const boxconstraints = BoxConstrainedSpace function __init__() @require UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" begin diff --git a/src/operators/crossover.jl b/src/operators/crossover.jl deleted file mode 100644 index 6d3f5377..00000000 --- a/src/operators/crossover.jl +++ /dev/null @@ -1,150 +0,0 @@ -""" - UniformCrossover(;p = 0.5) - -Uniform crossover aka Binomial crossover suitable for binary representation. -""" -struct UniformCrossover - p::Float64 - UniformCrossover(;p=0.5) = new(p) -end - -function crossover(population, parameters::UniformCrossover) - n = length(population) ÷ 2 - offspring_A = positions(population[1:n]) - offspring_B = positions(population[n+1:2n]) - mask = rand(size(offspring_A)...) .<= parameters.p - tmp = copy(offspring_A[mask]) - offspring_A[mask] = offspring_B[mask] - offspring_B[mask] = tmp - [offspring_A; offspring_B] -end - -""" - OrderCrossover() -Order crossover for representation where order is important. Suitable for permutation -representation. -""" -struct OrderCrossover end - -function crossover(population, parameters::OrderCrossover) - O = positions(population) - N, D = size(O) - s = rand(1:D, N) # slash points - for i = 1:2:N - PA = O[i, :] # parent A - PB = O[i+1, :] # parent B - O[i, s[i]+1:D] = setdiff(PB, PA[1:s[i]]); - O[i+1,s[i]+1:D] = setdiff(PA, PB[1:s[i]]); - end - O -end - -########################################################################## - -function gen_β(β, η, D, R) - α = 2.0 .- β .^ (- η - 1.0 ) - mask = R .<= 1.0 ./ α - s = 1.0 / (η + 1.0) - βq = [ mask[i] ? (R[i] * α[i])^s : (1.0 / (2.0 - R[i]*α[i]))^s for i in 1:D] - βq -end - -""" - SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9) - -Simulated binomial crossover for given two `Vectors{Real}`. -""" -function SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9) - xu = view(bounds, 2,:) - xl = view(bounds, 1,:) - D = length(vector1) - - do_crossover = ones(Bool, D) - do_crossover[rand(D) .> p_variable] .= false - do_crossover[ abs.( vector2 - vector1 ) .<= eps() ] .= false - - y1 = min.( vector1, vector2 ) - y2 = max.( vector1, vector2 ) - Δ = max.(eps(), y2 - y1) - - - R = rand(D) - - β = @. 1.0 + (2.0 * (y1 - xl) / Δ) - βq = gen_β(β, η, D, R) - c1 = @. 0.5*(y1 + y2 - βq*Δ) - - β = @. 1.0 + (2.0 * (y1 - xl) / Δ) - βq = gen_β(β, η, D, R) - c2 = @. 0.5*(y1 + y2 + βq*Δ) - - # swap - mask = rand(Bool, D) - cc = copy(c1) - c1[mask] = c2[mask] - c2[mask] = cc[mask] - - cc1 = copy(vector1) - cc1[do_crossover] = _to_int_if_necessary(eltype(cc1), c1[do_crossover] ) - cc2 = copy(vector2) - cc2[do_crossover] = _to_int_if_necessary(eltype(cc2), c2[do_crossover] ) - - - reset_to_violated_bounds!(cc1, bounds) - reset_to_violated_bounds!(cc2, bounds) - - return cc1, cc2 -end - -""" - SBX(;η, p, bounds) - -Simulated Binomial Crossover. -""" -mutable struct SBX - η::Float64 - p::Float64 - bounds::Matrix{Float64} - SBX(;η = 15, p = 0.9, bounds = zeros(0,0)) = new(η, p, bounds) -end - -function crossover(population, parameters::SBX) - isempty(population) && return zeros(0,0) - Q = positions(population) - bounds = parameters.bounds - for i in 1:2:length(population)-1 - p1 = get_position(population[i]) - p2 = get_position(population[i+1]) - c1, c2 = SBX_crossover(p1, p2, bounds, parameters.η, parameters.p) - Q[i,:] = c1 - Q[i+1,:] = c2 - end - Q -end - - - -""" - DE_crossover(x, u, CR) - -Binomial crossover between x and u for Differential Evolution with probability CR, i.e., -`v[j] = u[j]` if `rand() < CR`, otherwise `v[j] = x[j]`. Return `v`. -""" -function DE_crossover(x, u, CR) - D = length(x) - # binomial crossover - v = zeros(D) - j_rand = rand(1:D) - - # binomial crossover - for j = 1:D - if rand() < CR || j == j_rand - v[j] = u[j] - else - v[j] = x[j] - end - end - - return v -end - diff --git a/src/operators/crossover/order.jl b/src/operators/crossover/order.jl new file mode 100644 index 00000000..8f7412d4 --- /dev/null +++ b/src/operators/crossover/order.jl @@ -0,0 +1,23 @@ +""" + OrderCrossover() +Order crossover for representation where order is important. Suitable for permutation +representation. +""" +struct OrderCrossover + rng +end + +OrderCrossover(;rng = default_rng_mh()) = OrderCrossover(rng) + +function crossover(population, parameters::OrderCrossover) + O = positions(population) + N, D = size(O) + s = rand(parameters.rng, 1:D, N) # slash points + for i = 1:2:N + PA = O[i, :] # parent A + PB = O[i+1, :] # parent B + O[i, s[i]+1:D] = setdiff(PB, PA[1:s[i]]); + O[i+1,s[i]+1:D] = setdiff(PA, PB[1:s[i]]); + end + O +end diff --git a/src/operators/crossover/sbx.jl b/src/operators/crossover/sbx.jl new file mode 100644 index 00000000..36c560f7 --- /dev/null +++ b/src/operators/crossover/sbx.jl @@ -0,0 +1,105 @@ +""" + SBX(;η, p, bounds) + +Simulated Binomial Crossover. +""" +mutable struct SBX + η::Float64 + p::Float64 + bounds + rng +end + +SBX(;η=15, p=0.9, bounds=zeros(0,0), rng = default_rng_mh()) = SBX(η, p, bounds, rng) + + +function crossover(population, parameters::SBX) + isempty(population) && return zeros(0,0) + Q = positions(population) + bounds = parameters.bounds + rng = parameters.rng + for i in 1:2:length(population)-1 + p1 = get_position(population[i]) + p2 = get_position(population[i+1]) + c1, c2 = SBX_crossover(p1, p2, bounds, parameters.η, parameters.p, rng) + Q[i,:] = c1 + Q[i+1,:] = c2 + end + Q +end + +function gen_β(β, η, D, R) + α = 2 .- β .^ (- η - 1 ) + mask = R .<= 1 ./ α + s = 1 / (η + 1) + βq = [ mask[i] ? (R[i] * α[i])^s : (1 / (2- R[i]*α[i]))^s for i in 1:D] + βq +end + +""" + SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9) + +Simulated binomial crossover for given two `Vectors{Real}`. +""" +function SBX_crossover( + vector1, + vector2, + bounds::BoxConstrainedSpace, + η=15, + p_variable = 0.9, + rng = default_rng_mh() + ) + + xu = bounds.ub + xl = bounds.lb + D = getdim(bounds) + + do_crossover = ones(Bool, D) + do_crossover[rand(rng, D) .> p_variable] .= false + do_crossover[ abs.( vector2 - vector1 ) .<= eps() ] .= false + + y1 = min.( vector1, vector2 ) + y2 = max.( vector1, vector2 ) + Δ = max.(eps(), y2 - y1) + + + R = rand(rng, D) + + β = @. 1 + (2* (y1 - xl) / Δ) + βq = gen_β(β, η, D, R) + c1 = @. 0.5*(y1 + y2 - βq*Δ) + + β = @. 1 + (2* (y1 - xl) / Δ) + βq = gen_β(β, η, D, R) + c2 = @. 0.5*(y1 + y2 + βq*Δ) + + # swap + mask = rand(rng, Bool, D) + cc = copy(c1) + c1[mask] = c2[mask] + c2[mask] = cc[mask] + + cc1 = copy(vector1) + cc1[do_crossover] = _to_int_if_necessary(eltype(cc1), c1[do_crossover] ) + cc2 = copy(vector2) + cc2[do_crossover] = _to_int_if_necessary(eltype(cc2), c2[do_crossover] ) + + + reset_to_violated_bounds!(cc1, bounds) + reset_to_violated_bounds!(cc2, bounds) + + return cc1, cc2 +end + +function SBX_crossover( + vector1, + vector2, + bounds::AbstractMatrix, + η=15, + p_variable = 0.9, + rng = default_rng_mh() + ) + b = BoxConstrainedSpace(lb = bounds[1,:], ub = bounds[2,:]) + SBX_crossover(vector1, vector2, b, η, p_variable, rng) +end + diff --git a/src/operators/crossover/uniform.jl b/src/operators/crossover/uniform.jl new file mode 100644 index 00000000..5e95fa8c --- /dev/null +++ b/src/operators/crossover/uniform.jl @@ -0,0 +1,50 @@ +""" + UniformCrossover(;p = 0.5) + +Uniform crossover a.k.a. Binomial crossover. Suitable for binary representation. +""" +struct UniformCrossover + p::Float64 + rng +end + +UniformCrossover(;p=0.5, rng = default_rng_mh()) = UniformCrossover(p, rng) + +function crossover(population, parameters::UniformCrossover) + n = length(population) ÷ 2 + rng = parameters.rng + offspring_A = positions(population[1:n]) + offspring_B = positions(population[n+1:2n]) + mask = rand(rng, size(offspring_A)...) .<= parameters.p + # swapping genes + tmp = copy(offspring_A[mask]) + offspring_A[mask] = offspring_B[mask] + offspring_B[mask] = tmp + [offspring_A; offspring_B] +end + + +""" + DE_crossover(x, u, CR) + +Binomial crossover between x and u for Differential Evolution with probability CR, i.e., +`v[j] = u[j]` if `rand() < CR`, otherwise `v[j] = x[j]`. Return `v`. +""" +function DE_crossover(x, u, CR) + D = length(x) + # binomial crossover + v = zeros(D) + j_rand = rand(1:D) + + # binomial crossover + for j = 1:D + if rand() < CR || j == j_rand + v[j] = u[j] + else + v[j] = x[j] + end + end + + return v +end + diff --git a/src/operators/environmental_selection.jl b/src/operators/environmental_selection.jl index 2b855db3..b4d2ec21 100644 --- a/src/operators/environmental_selection.jl +++ b/src/operators/environmental_selection.jl @@ -29,3 +29,38 @@ function environmental_selection!( sort!(population, lt=is_better) deleteat!(population, N+1:length(population)) end + + +""" + RankAndCrowding() + +Perform `environmental_selection!` based non-dominated ranking and crowding +distance. +""" +struct RankAndCrowding end + +function environmental_selection!( + population, + offsprings, + parameters::RankAndCrowding; + is_better = is_better + ) + + N = length(population) + append!(population, offsprings) + fast_non_dominated_sort!(population) + + let f::Int = 1 + ind = 0 + indnext = findlast(x -> x.rank == f, population) + while !isnothing(indnext) && 0 < indnext <= N + ind = indnext + f += 1 + indnext = findlast(x -> x.rank == f, population) + end + isnothing(indnext) && (indnext = length(population)) + update_crowding_distance!(view(population, ind+1:indnext)) + sort!(view(population, ind+1:indnext), by = x -> x.crowding, rev = true, alg = PartialQuickSort(N-ind)) + end + deleteat!(population, N+1:length(population)) +end diff --git a/src/operators/initializer.jl b/src/operators/initializer.jl index 4cce5889..8c57d764 100644 --- a/src/operators/initializer.jl +++ b/src/operators/initializer.jl @@ -1,141 +1,5 @@ abstract type AbstractInitializer end -struct RandomInBounds <: AbstractInitializer - N::Int -end - -""" - RandomInBounds - -Initialize `N` solutions with random values in bounds. Suitable for -integer and real coded problems. -""" -RandomInBounds(;N = 0) = RandomInBounds(N) - -""" - RandomBinary(;N) - -Create random binary individuals. -""" -struct RandomBinary <: AbstractInitializer - N::Int - RandomInBounds(;N = 0) = new(N) -end - - - - -""" - RandomPermutation(;N) - -Create individuals in random permutations. -""" -struct RandomPermutation <: AbstractInitializer - N::Int - RandomPermutation(;N = 0) = new(N) -end - - -""" - LatinHypercubeSampling(nsamples, dim; iterations) - -Create `N` solutions within a Latin Hypercube sample in bounds with dim. - -### Example - -```julia-repl -julia> sample(LatinHypercubeSampling(10,2)) -10×2 Matrix{Float64}: - 0.0705631 0.795046 - 0.7127 0.0443734 - 0.118018 0.114347 - 0.48839 0.903396 - 0.342403 0.470998 - 0.606461 0.275709 - 0.880482 0.89515 - 0.206142 0.321041 - 0.963978 0.527518 - 0.525742 0.600209 - -julia> sample(LatinHypercubeSampling(10,2), [-10 -10;10 10.0]) -10×2 Matrix{Float64}: - -7.81644 -2.34461 - 0.505902 0.749366 - 3.90738 -8.57816 - -2.05837 9.803 - 5.62434 6.82463 - -9.34437 2.72363 - 6.43987 -1.74596 - -1.3162 -4.50273 - 9.45114 -7.13632 - -4.71696 5.0381 -``` -""" -struct LatinHypercubeSampling <: AbstractInitializer - N::Int - dim::Int - iterations::Int - LatinHypercubeSampling(nsamples, dim;iterations=25) = new(nsamples,dim,iterations) -end - -""" - Grid(npartitions, dim) - -Parameters to generate a grid with `npartitions` in a space with `dim` dimensions. - -### Example - -```julia-repl -julia> sample(Grid(5,2)) -25×2 Matrix{Float64}: - 0.0 0.0 - 0.25 0.0 - 0.5 0.0 - 0.75 0.0 - ⋮ - 0.5 1.0 - 0.75 1.0 - 1.0 1.0 - -julia> sample(Grid(5,2), [-1 -1; 1 1.]) -25×2 Matrix{Float64}: - -1.0 -1.0 - -0.5 -1.0 - 0.0 -1.0 - 0.5 -1.0 - ⋮ - 0.0 1.0 - 0.5 1.0 - 1.0 1.0 -``` - -Note that the sample is with size `npartitions^(dim)`. -""" -struct Grid <: AbstractInitializer - npartitions::Int - dim::Int -end - -function gen_initial_state(problem,parameters::RandomPermutation,information,options) - - D = size(problem.bounds, 2) - X = zeros(Int, parameters.N, D) - N = parameters.N - for i in 1:parameters.N - X[i,:] = shuffle(1:D) - end - - if problem.parallel_evaluation - population = create_solutions(X, problem; ε=options.h_tol) - else - population = [ create_solution(X[i,:], problem; ε=options.h_tol) for i in 1:N] - end - - best_solution = get_best(population) - - State(best_solution, population; f_calls = length(population), iteration=1) -end - function initialize!( status, parameters::AbstractInitializer, @@ -149,5 +13,3 @@ function initialize!( gen_initial_state(problem,parameters,information,options, status) end -include("sample.jl") - diff --git a/src/operators/mutation.jl b/src/operators/mutation.jl deleted file mode 100644 index f86a559a..00000000 --- a/src/operators/mutation.jl +++ /dev/null @@ -1,181 +0,0 @@ -""" - BitFlipMutation(;p = 1e-2) - -Flip each bit with probability `p`. -""" -struct BitFlipMutation - p::Float64 - BitFlipMutation(;p = 1e-2) = new(p) -end - -function mutation!(Q, parameters::BitFlipMutation) - mask = rand(size(Q)...) .< parameters.p - Q[mask] = .!Q[mask] - Q -end - -""" - SlightMutation - -> Fogel, D. B. (1988). An evolutionary approach to the traveling salesman problem. -> Biological Cybernetics, 60(2), 139-144. -""" -struct SlightMutation end - -function mutation!(Q, parameters::SlightMutation) - N, D = size(Q) - k = rand(1:D, N) - s = rand(1:D, N) - for i in 1:N - s[i] == k[i] && continue - if s[i] < k[i] - c = vcat(1:s[i]-1, k[i], s[i]:k[i]-1, k[i]+1:D) - else s[i] > k[i] - c = vcat(1:k[i]-1, k[i]+1:s[i]-1, k[i], s[i]:D) - end - Q[i,:] = Q[i, c] - end - Q -end - - -###################### - -""" - polynomial_mutation!(vector, bounds, η=20, prob = 1 / length(vector)) - -Polynomial Mutation applied to a vector of real numbers. -""" -function polynomial_mutation!(vector, bounds, η=20, prob = 1 / length(vector)) - do_mutation = rand(length(vector)) .< prob - - xu = view(bounds, 2,do_mutation) - xl = view(bounds, 1,do_mutation) - x = view(vector, do_mutation) - - δ1 = (x - xl) ./ (xu - xl) - δ2 = (xu - x) ./ (xu - xl) - - D = length(xu) - R = rand(D) - mask = R .< 0.5 - s = η+1.0 - mut_pow = 1.0 / (η + 1.0) - δq = [ mask[i] ? - ^(2.0R[i] + (1. - 2.0R[i]) * ^(1.0 - δ1[i], s), mut_pow) - 1.0 : - 1.0 - (2.0 * (1.0 - R[i]) + 2.0 * (R[i] - 0.5) * ^(1.0 - δ2[i], s))^mut_pow - for i in 1:D - ] - - vector[do_mutation] = _to_int_if_necessary(eltype(vector), x + δq .* ( xu - xl)) - vector -end - -""" - PolynomialMutation(;η, p, bounds) - -Polynomial mutation. -""" -struct PolynomialMutation - η::Float64 - p::Float64 - bounds::Matrix{Float64} - PolynomialMutation(;η = 15.0, bounds = zeros(0,0), p=1/size(bounds,2)) = begin - new(η, isfinite(p) ? p : 1e-2, bounds) - end -end - -function mutation!(Q, parameters::PolynomialMutation) - for i in 1:size(Q,1) - polynomial_mutation!(view(Q, i,:), parameters.bounds, parameters.η, parameters.p) - end - Q -end - -""" - DE_mutation(population, F = 1.0, strategy = :rand1) - -Generate a `Vector` computed from population used in Differential Evolution. -Parameters: `F` is the stepsize, `strategy` can be one the following `:best1`, `:rand2`, -`:randToBest1` or `:best2`. -""" -function DE_mutation(population, - F::Float64 = 1.0, - strategy::Symbol=:rand1, - best_ind=0 - ) - - N = length(population) - - i = rand(1:N) - # select participats - r1 = rand(1:N) - while r1 == i - r1 = rand(1:N) - end - - r2 = rand(1:N) - while r2 == i || r1 == r2 - r2 = rand(1:N) - end - - - x = get_position(population[i] ) - a = get_position(population[r1]) - b = get_position(population[r2]) - - # strategy is selected here - if strategy == :rand1 - # DE/rand/1 - return x + F * (a - b) - end - - r3 = rand(1:N) - while r3 == i || r3 == r1 || r3 == r2 - r3 = rand(1:N) - end - - c = get_position(population[r3]) - - if strategy == :rand2 - # DE/rand/2 - - r4 = rand(1:N) - while r4 == i || r4 == r1 || r4 == r2 || r4 == r3 - r4 = rand(1:N) - end - - r5 = rand(1:N) - while r5 == i || r5 == r1 || r5 == r2 || r5 == r3 || r5 == r4 - r5 = rand(1:N) - end - - d = get_position(population[r4]) - ee = get_position(population[r5]) - - return ee + F * (a - b + c - d) - end - - best_ind = best_ind == 0 ? argbest(population) : best_ind - x_best = get_position(population[best_ind]) - if strategy == :best1 - # DE/best/1 - u = x_best + F * (b - c) - elseif strategy == :randToBest1 - # DE/rand-to-best/1 - u = x + F * (x_best - x + a - b) - elseif strategy == :best2 - # DE/best/2 - r4 = rand(1:N) - while r4 == i || r4 == r1 || r4 == r2 || r4 == r3 || r4 == best_ind - r4 = rand(1:N) - end - d = population[r4].x - u = x_best + F * (a - b + c - d) - else - error("Unknown strategy " * string(strategy)) - end - - return u -end - diff --git a/src/operators/mutation/bitflip.jl b/src/operators/mutation/bitflip.jl new file mode 100644 index 00000000..2cb6e2ad --- /dev/null +++ b/src/operators/mutation/bitflip.jl @@ -0,0 +1,17 @@ +""" + BitFlipMutation(;p = 1e-2) + +Flip each bit with probability `p`. +""" +struct BitFlipMutation + p::Float64 + rng +end + +BitFlipMutation(;p = 1e-2, rng = default_rng_mh()) = BitFlipMutation(p, rng) + +function mutation!(Q::AbstractMatrix{Bool}, parameters::BitFlipMutation) + mask = rand(parameters.rng, size(Q)...) .< parameters.p + Q[mask] = .!Q[mask] + Q +end diff --git a/src/operators/mutation/mutation.jl b/src/operators/mutation/mutation.jl new file mode 100644 index 00000000..e9b6c502 --- /dev/null +++ b/src/operators/mutation/mutation.jl @@ -0,0 +1,87 @@ +""" + DE_mutation(population, F = 1.0, strategy = :rand1) + +Generate a `Vector` computed from population used in Differential Evolution. +Parameters: `F` is the stepsize, `strategy` can be one the following `:best1`, `:rand2`, +`:randToBest1` or `:best2`. +""" +function DE_mutation(population, + F::Float64 = 1.0, + strategy::Symbol=:rand1, + best_ind=0 + ) + + N = length(population) + + i = rand(1:N) + # select participats + r1 = rand(1:N) + while r1 == i + r1 = rand(1:N) + end + + r2 = rand(1:N) + while r2 == i || r1 == r2 + r2 = rand(1:N) + end + + + x = get_position(population[i] ) + a = get_position(population[r1]) + b = get_position(population[r2]) + + # strategy is selected here + if strategy == :rand1 + # DE/rand/1 + return x + F * (a - b) + end + + r3 = rand(1:N) + while r3 == i || r3 == r1 || r3 == r2 + r3 = rand(1:N) + end + + c = get_position(population[r3]) + + if strategy == :rand2 + # DE/rand/2 + + r4 = rand(1:N) + while r4 == i || r4 == r1 || r4 == r2 || r4 == r3 + r4 = rand(1:N) + end + + r5 = rand(1:N) + while r5 == i || r5 == r1 || r5 == r2 || r5 == r3 || r5 == r4 + r5 = rand(1:N) + end + + d = get_position(population[r4]) + ee = get_position(population[r5]) + + return ee + F * (a - b + c - d) + end + + best_ind = best_ind == 0 ? argbest(population) : best_ind + x_best = get_position(population[best_ind]) + if strategy == :best1 + # DE/best/1 + u = x_best + F * (b - c) + elseif strategy == :randToBest1 + # DE/rand-to-best/1 + u = x + F * (x_best - x + a - b) + elseif strategy == :best2 + # DE/best/2 + r4 = rand(1:N) + while r4 == i || r4 == r1 || r4 == r2 || r4 == r3 || r4 == best_ind + r4 = rand(1:N) + end + d = population[r4].x + u = x_best + F * (a - b + c - d) + else + error("Unknown strategy " * string(strategy)) + end + + return u +end + diff --git a/src/operators/mutation/polynomial.jl b/src/operators/mutation/polynomial.jl new file mode 100644 index 00000000..d15fc5a3 --- /dev/null +++ b/src/operators/mutation/polynomial.jl @@ -0,0 +1,81 @@ +""" + PolynomialMutation(;η, p, bounds) + +Polynomial mutation. +""" +struct PolynomialMutation + η::Float64 + p::Float64 + bounds + rng +end + +function PolynomialMutation(; + η = 15.0, + bounds = zeros(0,0), + p=1/getdim(bounds), + rng = default_rng_mh() + ) + PolynomialMutation(η, isfinite(p) ? p : 1e-2, bounds, rng) +end + +function mutation!(Q, parameters::PolynomialMutation) + bounds = parameters.bounds + η = parameters.η + p = parameters.p + rng = parameters.rng + for i in 1:size(Q,1) + polynomial_mutation!(view(Q, i,:), bounds, η, p, rng) + end + Q +end + + +""" + polynomial_mutation!(vector, bounds, η=20, prob = 1 / length(vector)) + +Polynomial Mutation applied to a vector of real numbers. +""" +function polynomial_mutation!( + vector, + bounds::BoxConstrainedSpace, + η=20, + prob = 1 / length(vector), + rng = default_rng_mh() + ) + + do_mutation = rand(rng, length(vector)) .< prob + + xu = view(bounds.ub, do_mutation) + xl = view(bounds.lb, do_mutation) + x = view(vector, do_mutation) + + δ1 = (x - xl) ./ (xu - xl) + δ2 = (xu - x) ./ (xu - xl) + + D = length(xu) + R = rand(rng, D) + mask = R .< 0.5 + s = η+1 + mut_pow = 1 / (η + 1) + δq = [ mask[i] ? + ^(2R[i] + (1. - 2R[i]) * ^(1 - δ1[i], s), mut_pow) - 1 : + 1 - (2 * (1 - R[i]) + 2 * (R[i] - 0.5) * ^(1 - δ2[i], s))^mut_pow + for i in 1:D + ] + + vector[do_mutation] = _to_int_if_necessary(eltype(vector), x + δq .* ( xu - xl)) + vector +end + +function polynomial_mutation!(vector, + bounds::AbstractMatrix, + η=20, + prob = 1 / length(vector), + rng = default_rng_mh() + ) + + b = BoxConstrainedSpace(lb = bounds[1,:], ub = bounds[2,:]) + polynomial_mutation!(vector, b, η, prob, rng) +end + diff --git a/src/operators/mutation/slight.jl b/src/operators/mutation/slight.jl new file mode 100644 index 00000000..4e56ffb7 --- /dev/null +++ b/src/operators/mutation/slight.jl @@ -0,0 +1,28 @@ +""" + SlightMutation + +> Fogel, D. B. (1988). An evolutionary approach to the traveling salesman problem. +> Biological Cybernetics, 60(2), 139-144. +""" +struct SlightMutation + rng +end + +SlightMutation(;rng = default_rng_mh()) = SlightMutation(rng) + +function mutation!(Q, parameters::SlightMutation) + N, D = size(Q) + k = rand(parameters.rng, 1:D, N) + s = rand(parameters.rng, 1:D, N) + for i in 1:N + s[i] == k[i] && continue + if s[i] < k[i] + c = vcat(1:s[i]-1, k[i], s[i]:k[i]-1, k[i]+1:D) + else s[i] > k[i] + c = vcat(1:k[i]-1, k[i]+1:s[i]-1, k[i], s[i]:D) + end + Q[i,:] = Q[i, c] + end + Q +end + diff --git a/src/operators/operators.jl b/src/operators/operators.jl index 4beaa36b..f8269702 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -1,6 +1,21 @@ include("utils.jl") + include("initializer.jl") -include("selection.jl") -include("crossover.jl") -include("mutation.jl") +include("sampling/random.jl") +include("sampling/grid.jl") +include("sampling/latinhypercube.jl") +include("sampling/searchspaces.jl") + +include("selection/selection.jl") + +# crossover +include("crossover/order.jl") +include("crossover/sbx.jl") +include("crossover/uniform.jl") + +include("mutation/bitflip.jl") +include("mutation/polynomial.jl") +include("mutation/slight.jl") +include("mutation/mutation.jl") + include("environmental_selection.jl") diff --git a/src/operators/sampling/grid.jl b/src/operators/sampling/grid.jl new file mode 100644 index 00000000..31e98887 --- /dev/null +++ b/src/operators/sampling/grid.jl @@ -0,0 +1,51 @@ +""" + Grid(npartitions, dim) + +Parameters to generate a grid with `npartitions` in a space with `dim` dimensions. + +### Example + +```julia-repl +julia> sample(Grid(5,2)) +25×2 Matrix{Float64}: + 0.0 0.0 + 0.25 0.0 + 0.5 0.0 + 0.75 0.0 + ⋮ + 0.5 1.0 + 0.75 1.0 + 1.0 1.0 + +julia> sample(Grid(5,2), [-1 -1; 1 1.]) +25×2 Matrix{Float64}: + -1.0 -1.0 + -0.5 -1.0 + 0.0 -1.0 + 0.5 -1.0 + ⋮ + 0.0 1.0 + 0.5 1.0 + 1.0 1.0 +``` + +Note that the sample is with size `npartitions^(dim)`. +""" +struct Grid <: AbstractInitializer + npartitions::Int + dim::Int +end + +function sample(method::Grid, bounds = zeros(0,0)) + v = range(0, 1, length = method.npartitions) + vals = Iterators.product((v for _ in 1:method.dim)...) + + # fill values in the grid + _X = reshape([x for val in vals for x in val], method.dim, length(vals)) + X = Array(_X') + + isempty(bounds) && (return X) + _scale_sample(X, bounds) +end + + diff --git a/src/operators/sample.jl b/src/operators/sampling/latinhypercube.jl similarity index 64% rename from src/operators/sample.jl rename to src/operators/sampling/latinhypercube.jl index b65ebd4e..149b54b8 100644 --- a/src/operators/sample.jl +++ b/src/operators/sampling/latinhypercube.jl @@ -1,17 +1,43 @@ -function _lhs(nsamples, dim) - # initial sample - X = reshape([v for i in 1:dim for v in shuffle(1.0:nsamples)], nsamples, dim) - # smooth and normalize sample - (X - rand(nsamples, dim)) / nsamples -end +""" + LatinHypercubeSampling(nsamples, dim; iterations) -_score_lhs(M) = minimum(pairwise_distances(M)) -function _scale_sample(X, bounds) - a = view(bounds, 1,:)' - b = view(bounds, 2,:)' - @assert length(a) == size(X,2) - # scale sample - a .+ (b - a) .* X +Create `N` solutions within a Latin Hypercube sample in bounds with dim. + +### Example + +```julia-repl +julia> sample(LatinHypercubeSampling(10,2)) +10×2 Matrix{Float64}: + 0.0705631 0.795046 + 0.7127 0.0443734 + 0.118018 0.114347 + 0.48839 0.903396 + 0.342403 0.470998 + 0.606461 0.275709 + 0.880482 0.89515 + 0.206142 0.321041 + 0.963978 0.527518 + 0.525742 0.600209 + +julia> sample(LatinHypercubeSampling(10,2), [-10 -10;10 10.0]) +10×2 Matrix{Float64}: + -7.81644 -2.34461 + 0.505902 0.749366 + 3.90738 -8.57816 + -2.05837 9.803 + 5.62434 6.82463 + -9.34437 2.72363 + 6.43987 -1.74596 + -1.3162 -4.50273 + 9.45114 -7.13632 + -4.71696 5.0381 +``` +""" +struct LatinHypercubeSampling <: AbstractInitializer + N::Int + dim::Int + iterations::Int + LatinHypercubeSampling(nsamples, dim;iterations=25) = new(nsamples,dim,iterations) end """ @@ -66,17 +92,19 @@ function sample(method::LatinHypercubeSampling, bounds = zeros(0,0)) _scale_sample(X, bounds) end -function sample(method::Grid, bounds = zeros(0,0)) - v = range(0, 1, length = method.npartitions) - vals = Iterators.product((v for _ in 1:method.dim)...) - - # fill values in the grid - _X = reshape([x for val in vals for x in val], method.dim, length(vals)) - X = Array(_X') - - isempty(bounds) && (return X) - _scale_sample(X, bounds) +function _lhs(nsamples, dim) + # initial sample + X = reshape([v for i in 1:dim for v in shuffle(1.0:nsamples)], nsamples, dim) + # smooth and normalize sample + (X - rand(nsamples, dim)) / nsamples end +_score_lhs(M) = minimum(pairwise_distances(M)) +function _scale_sample(X, bounds) + a = view(bounds, 1,:)' + b = view(bounds, 2,:)' + @assert length(a) == size(X,2) + # scale sample + a .+ (b - a) .* X +end -sample(method::RandomInBounds, bounds) = _scale_sample(rand(method.N, size(bounds,2)), bounds) diff --git a/src/operators/sampling/random.jl b/src/operators/sampling/random.jl new file mode 100644 index 00000000..812da91d --- /dev/null +++ b/src/operators/sampling/random.jl @@ -0,0 +1,61 @@ +struct RandomInBounds <: AbstractInitializer + N::Int +end + +""" + RandomInBounds + +Initialize `N` solutions with random values in bounds. Suitable for +integer and real coded problems. +""" +RandomInBounds(;N = 0) = RandomInBounds(N) + + +sample(method::RandomInBounds, bounds) = _scale_sample(rand(method.N, size(bounds,2)), bounds) + +""" + RandomBinary(;N) + +Create random binary individuals. +""" +struct RandomBinary <: AbstractInitializer + N::Int + RandomInBounds(;N = 0) = new(N) +end + + + + +""" + RandomPermutation(;N) + +Create individuals in random permutations. +""" +struct RandomPermutation <: AbstractInitializer + N::Int + RandomPermutation(;N = 0) = new(N) +end + +function gen_initial_state(problem,parameters::RandomPermutation,information,options) + + D = getdim(problem) + N = parameters.N + + if problem.search_space isa BoxConstrainedSpace + search_space = PermutationSpace(getdim(problem.search_space)) + else + search_space = problem.search_space + end + + X = sample(RandomSampler(search_space; options.rng), N) + + if problem.parallel_evaluation + population = create_solutions(X, problem; ε=options.h_tol) + else + population = [ create_solution(X[i,:], problem; ε=options.h_tol) for i in 1:N] + end + + best_solution = get_best(population) + + State(best_solution, population; f_calls = length(population), iteration=1) +end diff --git a/src/operators/sampling/searchspaces.jl b/src/operators/sampling/searchspaces.jl new file mode 100644 index 00000000..01493262 --- /dev/null +++ b/src/operators/sampling/searchspaces.jl @@ -0,0 +1,41 @@ +function sample( + sampler::Sampler{M, S}, + n = 100 + ) where {M <: RandomSampler, S <: AtomicSearchSpace} + + X = rand(sampler.method.rng, sampler.searchspace, n) + # convert to matrix + [X[i][j] for i in eachindex(X), j in eachindex(first(X))] +end + + +function sample( + sampler::Sampler{M, S}, args... + ) where {M <: GridSampler, S <: AtomicSearchSpace} + + X = collect(sampler) + # convert to matrix + [X[i][j] for i in eachindex(X), j in eachindex(first(X))] +end + +function sample( + ::Type{T}, + search_space::AtomicSearchSpace, + n = 100, + ) where T <: AbstractSampler + + sample(T(search_space), n) +end + +function sample(search_space::AtomicSearchSpace, n = 100) + + if cardinality(search_space) <= n + sampler = GridSampler(search_space) + else + sampler = RandomSampler(search_space) + end + + sample(sampler, n) +end + + diff --git a/src/operators/selection.jl b/src/operators/selection/selection.jl similarity index 100% rename from src/operators/selection.jl rename to src/operators/selection/selection.jl diff --git a/src/optimize.jl b/src/optimize.jl deleted file mode 100644 index 20fdc252..00000000 --- a/src/optimize.jl +++ /dev/null @@ -1,227 +0,0 @@ -""" - optimize( - f::Function, # objective function - bounds::AbstractMatrix, - method::AbstractAlgorithm = ECA(); - logger::Function = (status) -> nothing, - ) - -Minimize a n-dimensional function `f` with domain `bounds` (2×n matrix) using `method = ECA()` by default. - -# Example -Minimize f(x) = Σx² where x ∈ [-10, 10]³. - -Solution: - -```jldoctest -julia> f(x) = sum(x.^2) -f (generic function with 1 method) - -julia> bounds = [ -10.0 -10 -10; # lower bounds - 10.0 10 10 ] # upper bounds -2×3 Array{Float64,2}: - -10.0 -10.0 -10.0 - 10.0 10.0 10.0 - -julia> result = optimize(f, bounds) -+=========== RESULT ==========+ - iteration: 1429 - minimum: 2.5354499999999998e-222 - minimizer: [-1.5135301653303966e-111, 3.8688354844737692e-112, 3.082095708730726e-112] - f calls: 29989 - total time: 0.1543 s -+============================+ -``` -""" -function optimize( - f::Function, # objective function - bounds::AbstractMatrix, - method::AbstractAlgorithm = ECA(); - logger::Function = (status) -> nothing, - ) - - - - ##################################### - # common methods - ##################################### - # status = method.status - information = method.information - options = method.options - parameters = method.parameters - ################################### - - problem = Problem(f, Array(bounds); parallel_evaluation=options.parallel_evaluation) - seed!(options.seed) - - ################################### - - start_time = time() - - status = method.status - options.debug && @info("Initializing population...") - status = initialize!(status,parameters, problem, information, options) - status.start_time = start_time - method.status = status - - show_status(status, parameters, options) - - status.iteration = 1 - - - convergence = State{typeof(status.best_sol)}[] - - - - ################################### - # store convergence - ################################### - if options.store_convergence - update_convergence!(convergence, status) - end - - options.debug && @info("Starting main loop...") - - logger(status) - - status.stop = status.stop || default_stop_check(status, information, options) - while !status.stop - status.iteration += 1 - - update_state!(status, parameters, problem, information, options) - - # store the number of fuction evaluations - status.f_calls = problem.f_calls - - show_status(status, parameters, options) - - if options.store_convergence - update_convergence!(convergence, status) - end - - status.overall_time = time() - status.start_time - logger(status) - - # common stop criteria - status.stop = status.stop || default_stop_check(status, information, options) - - # user defined stop criteria - stop_criteria!(status, parameters, problem, information, options) - - end - - status.overall_time = time() - status.start_time - - final_stage!( - status, - parameters, - problem, - information, - options - ) - - options.debug && @info "Termination reason: " * termination_status_message(status) - status.convergence = convergence - - return status - -end - - -function show_status(status, parameters, options) - !options.debug && (return) - status.final_time = time() - msg = "Current Status of " * string(typeof(parameters)) - @info msg - display(status) -end - -""" - optimize!(f, bounds, method;logger) - -Perform an iteration of `method`, and save the results in `method.status`. - -### Example - -```julia -f, bounds, _ = Metaheuristics.TestProblems.sphere(); -method = ECA() -while !Metaheuristics.should_stop(method) - optimize!(f, bounds, method) -end -result = Metaheuristics.get_result(method) -``` - -See also [`optimize`](@docs). -""" -function optimize!( - f::Function, # objective function - bounds::AbstractMatrix, - method::AbstractAlgorithm; - logger::Function = (status) -> nothing, - ) - - ##################################### - # common settings - ##################################### - information = method.information - options = method.options - parameters = method.parameters - problem = Problem(f, Array(bounds); parallel_evaluation=options.parallel_evaluation) - ################################### - - status = method.status - problem.f_calls = status.f_calls - - # initialization - if status.iteration == 0 || isnothing(status.best_sol) - options.debug && @info("Initializing population...") - seed!(options.seed) - start_time = time() - status = initialize!(status,parameters, problem, information, options) - status.start_time = start_time - method.status = status - status.final_time = time() - show_status(status, parameters, options) - return method - end - - convergence = State{typeof(status.best_sol)}[] - - status.iteration += 1 - - update_state!(status, parameters, problem, information, options) - # store the number of fuction evaluations - status.f_calls = problem.f_calls - status.overall_time = time() - status.start_time - status.final_time = time() - - if options.store_convergence - update_convergence!(convergence, status) - end - - show_status(status, parameters, options) - logger(status) - - # common stop criteria - status.stop = status.stop || default_stop_check(status, information, options) - # user defined stop criteria - stop_criteria!(status, parameters, problem, information, options) - - - if status.stop - options.debug && @info "Performing final stage due to stop criteria." - final_stage!( - status, - parameters, - problem, - information, - options - ) - end - - - options.debug && status.stop && @info "Should stop because: " * termination_status_message(status) - - return method -end diff --git a/src/optimize/after.jl b/src/optimize/after.jl new file mode 100644 index 00000000..8fdfb8ba --- /dev/null +++ b/src/optimize/after.jl @@ -0,0 +1,11 @@ +function _after_optimization!(status, parameters, problem, information, options, convergence) + + + t = time() + final_stage!(status, parameters, problem, information, options) + status.overall_time += time() - t + + options.debug && @info "Termination reason: " * termination_status_message(status) + # save convergence + status.convergence = convergence; +end diff --git a/src/optimize/before.jl b/src/optimize/before.jl new file mode 100644 index 00000000..72cc2766 --- /dev/null +++ b/src/optimize/before.jl @@ -0,0 +1,61 @@ +function _before_optimization!(f, search_space, method, logger) + # check whether optimizer is compatible with search_space + check_compat(search_space, method.parameters) + + # TODO: method = deepcopy(method) + information = method.information + options = method.options + parameters = method.parameters + problem = Problem(f, search_space; parallel_evaluation=options.parallel_evaluation) + seed!(options.seed) + + # initialization steps + status = method.status + options.debug && @info("Initializing population...") + start_time = time() + # initialize population (... and optimizer) + status = initialize!(status, parameters, problem, information, options) + + # save times + status.start_time = start_time + status.overall_time += time() - start_time + method.status = status + termination = method.termination + + # set termination based on convergence if necessary + if isempty(termination.checkall) && isempty(termination.checkany) + if fval(first(status.population)) isa Array + # multi-objective termination criterion + c = RobustConvergence(ftol=options.f_tol) + options.debug && @info("Set termination criteria: RobustConvergence") + else + # single-objective criteria (convergence based) + c = CheckConvergence(f_tol_abs = options.f_tol, + f_tol_rel = options.f_tol_rel, + x_tol = options.x_tol) + options.debug && @info("Set termination criteria: convergence indicators") + end + + push!(termination.checkany, c) + end + + + # show status if necessary + show_status(status, parameters, options) + + status.iteration = 1 + convergence = State{typeof(status.best_sol)}[] + options.store_convergence && update_convergence!(convergence, status) + + options.debug && @info("Starting main loop...") + logger(status) + + status.stop = status.stop || + # common stop criteria (budget based) + default_stop_check(status, information, options) || + # stop criteria given by user + stop_check(status, termination) + + status, parameters, problem, information, options, convergence +end + diff --git a/src/optimize/during.jl b/src/optimize/during.jl new file mode 100644 index 00000000..96008d07 --- /dev/null +++ b/src/optimize/during.jl @@ -0,0 +1,24 @@ +function _during_optimization!(status, parameters, problem, information, options, convergence, logger, termination) + status.iteration += 1 + t = time() + # perform one optimization step + update_state!(status, parameters, problem, information, options) + status.overall_time += time() - t + # store the number of function evaluations + status.f_calls = problem.f_calls + # show status if debug = true + show_status(status, parameters, options) + # store convergence if necessary + options.store_convergence && update_convergence!(convergence, status) + # user defined logger + logger(status) + status.stop = status.stop || + # common stop criteria (budget based) + default_stop_check(status, information, options) || + # stop criteria given by user + stop_check(status, termination) + + # TODO: check if following method should be deprecated in v4 + # stop criteria defined by optimizer + !status.stop && stop_criteria!(status, parameters, problem, information, options) +end diff --git a/src/optimize/optimize.jl b/src/optimize/optimize.jl new file mode 100644 index 00000000..e7c2bb31 --- /dev/null +++ b/src/optimize/optimize.jl @@ -0,0 +1,151 @@ +include("utils.jl") +include("before.jl") +include("during.jl") +include("after.jl") + + +""" + optimize( + f::Function, # objective function + search_space, + method::AbstractAlgorithm = ECA(); + logger::Function = (status) -> nothing, + ) + +Minimize a n-dimensional function `f` with domain `search_space` (2×n matrix) using `method = ECA()` by default. + +# Example +Minimize f(x) = Σx² where x ∈ [-10, 10]³. + +Solution: + +```jldoctest +julia> f(x) = sum(x.^2) +f (generic function with 1 method) + +julia> bounds = [ -10.0 -10 -10; # lower bounds + 10.0 10 10 ] # upper bounds +2×3 Array{Float64,2}: + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + +julia> result = optimize(f, bounds) ++=========== RESULT ==========+ + iteration: 1429 + minimum: 2.5354499999999998e-222 + minimizer: [-1.5135301653303966e-111, 3.8688354844737692e-112, 3.082095708730726e-112] + f calls: 29989 + total time: 0.1543 s ++============================+ +``` +""" +function optimize( + f::Function, # objective function + search_space, + method::AbstractAlgorithm = ECA(); + logger::Function = (status) -> nothing, + ) + + status, parameters, problem, information, options, convergence = _before_optimization!(f, search_space, method, logger) + + while !status.stop + _during_optimization!(status, + parameters, + problem, + information, + options, + convergence, + logger, + method.termination + ) + end + + _after_optimization!(status, parameters, problem, information, options, convergence) + + status +end + + +""" + optimize!(f, search_space, method;logger) + +Perform an iteration of `method`, and save the results in `method.status`. + +### Example + +```julia +f, bounds, _ = Metaheuristics.TestProblems.sphere(); +method = ECA() +while !Metaheuristics.should_stop(method) + optimize!(f, bounds, method) +end +result = Metaheuristics.get_result(method) +``` + +See also [`optimize`](@docs). +""" +function optimize!( + f::Function, # objective function + search_space, + method::AbstractAlgorithm; + logger::Function = (status) -> nothing, + ) + + status = method.status + # initialization + if status.iteration == 0 || isnothing(status.best_sol) + _before_optimization!(f, search_space, method, logger) + return method + end + + options = method.options + problem = Problem(f, search_space; + parallel_evaluation=options.parallel_evaluation) + problem.f_calls = status.f_calls + + # main optimization step + _during_optimization!( + status, + method.parameters, + problem, + method.information, + options, + empty(status.population), # convergence + logger, + method.termination + ) + + if status.stop + options.debug && @info "Performing final stage due to stop criteria." + _after_optimization!( + status, + method.parameters, + problem, + method.information, + options, + empty(status.population), # convergence + ) + end + + + options.debug && status.stop && @info "Should stop because: " * termination_status_message(status) + + return method +end + + +function optimize( + f::Function, + _search_space, + ::Type{T}; + logger::Function = (status) -> nothing, + kargs... + ) where T <: AbstractParameters + + search_space = _mat_to_bounds(_search_space) + # configure parameters depending on the search_space + algo = get_parameters(f, search_space, T) + set_user_parameters!(algo; kargs...) + # call optimize api + optimize(f, search_space, algo; logger) +end diff --git a/src/optimize/utils.jl b/src/optimize/utils.jl new file mode 100644 index 00000000..c11c7b6e --- /dev/null +++ b/src/optimize/utils.jl @@ -0,0 +1,85 @@ +""" + get_parameters(f, search_space, H) + +Get default parameters for metaheuristic `H` regarding f and the search_space. + +### Example + +```julia +ga = Metaheuristics.get_parameters(f, Permutations(10), GA) +``` +""" +get_parameters(f, search_space, ::Type{T}) where T <: AbstractParameters = T() + + +function show_status_oneline(status, parameters, options) + # sorry of the hard code :-) + !options.verbose && (return) + d = Any[ + "Iteration" => status.iteration, + "Num. Evals" => nfes(status), + ] + # show header + + t = status.iteration + e = nfes(status) + m = minimum(status) + if m isa Number + push!(d, " Minimum " => m) + else + n = length(get_non_dominated_solutions(status.population)) + s = sprint(print, "$n/$(length(status.population))") + push!(d, " NDS " => s) + # + end + # feas = count(is_feasible.(status.population)) ÷ length(status.population) + try + push!(d, " Mean CVio" => mean(sum_violations.(status.population))) + catch + end + + push!(d, " Time " => @sprintf("%.4f s", status.overall_time)) + + if m isa Number + v = stop_check(status, + CheckConvergence( + f_tol_abs = options.f_tol, + f_tol_rel = options.f_tol_rel, + x_tol =options.x_tol + ), + report = false + ) ? "Yes" : "No" + push!(d, "Converged " => @sprintf("%10s", v)) + end + + + if status.iteration <= 1 || status.iteration % 1000 == 0 + nm = first.(d) + lines = [fill('-', length(n) + 2) |> join for n in nm] + println("+", join(lines, "+"), "+") + println("| ", join(nm, " | "), " |") + println("+", join(lines, "+"), "+") + end + print("|") + + for v in last.(d) + if v isa Integer + @printf("% 10d | ", v) + elseif v isa AbstractString + @printf("% 10s | ", v) + elseif v isa AbstractFloat + @printf("%1.4e | ", v) + else + print(v, " | ") + end + end + println("") +end + +function show_status(status, parameters, options) + !options.debug && (return show_status_oneline(status, parameters, options)) + status.final_time = time() + msg = "Current Status of " * string(typeof(parameters)) + @info msg + display(status) +end diff --git a/src/parameters/parameters.jl b/src/parameters/parameters.jl new file mode 100644 index 00000000..88106a69 --- /dev/null +++ b/src/parameters/parameters.jl @@ -0,0 +1,21 @@ +function set_user_parameters!(algo::AbstractAlgorithm; kargs...) + parameters = algo.parameters + options = algo.options + # TODO information = algo.information + T = typeof(parameters) + for (field, value) in kargs + # set algorithm parameters + if ismutable(parameters) && hasfield(T, field) + setproperty!(parameters, field, value) + # set optional parameters + elseif ismutable(options) && hasfield(Options, field) + setproperty!(options, field, value) + # TODO elseif ismutable(information) && hasfield(Information, field) + # TODO setproperty!(information, field, value) + else + @warn "Parameter `$field = $value` never used." + end + end + algo +end + diff --git a/src/solutions/display.jl b/src/solutions/display.jl index 365281fc..06715153 100644 --- a/src/solutions/display.jl +++ b/src/solutions/display.jl @@ -3,18 +3,18 @@ import Base.show function print_vector(io, vector) n_items_in_vector = 5 if length(vector) < n_items_in_vector - show(io, vector) + show(IOContext(io, :compact => true), vector) return end @printf(io, "[") - for x in vector[1:2] - @printf(io, "%1.3e, ", x) + for x in vector[1:3] + @printf(io, "%g, ", x) end print(io, "…, ") for x in vector[end:end] - @printf(io, "%1.3e", x) + @printf(io, "%g", x) end @printf(io, "]") @@ -55,47 +55,50 @@ function Base.show(io::IO, solution::xFgh_indiv) end +show_pf(io, pf) = nothing function describe_result(io, status::State{T}) where T <: xFgh_solution - @printf(io, "%12s", "population:") - show(io, "text/plain", Array(status.population)) - # non-dominated pf = get_non_dominated_solutions(status.population) - println(io, "\nnon-dominated solution(s):") - show(io, "text/plain", pf) - print(io, "\n") + @printf(io," %-16s ", "Non-dominated:") + println(io, length(pf), " / ", length(status.population)) + show_pf(io, pf) end function describe_result(io, status::State) - @printf(io,"%12s %g\n", "minimum:", minimum(status)) - @printf(io,"%12s ", "minimizer:") - show(io, minimizer(status)) + @printf(io," %-16s %g\n", "Minimum:", minimum(status)) + @printf(io," %-16s ", "Minimizer:") + print_vector(io, minimizer(status)) println(io, "") end function Base.show(io::IO, status::State) - println(io, "+=========== RESULT ==========+") - @printf(io,"%12s %.0f\n", "iteration:", status.iteration) + title = "Optimization Result" + decor = join(repeat('=', length(title))) + printstyled(io, title, "\n", color=:blue, bold=true) + printstyled(io, decor, "\n", color=:gray) + + status.best_sol isa Nothing && return println(io, " Empty status.") + + @printf(io," %-16s %.0f\n", "Iteration:", status.iteration) describe_result(io, status) - @printf(io,"%12s %.0f\n", "f calls:", status.f_calls) + @printf(io," %-16s %.0f\n", "Function calls:", status.f_calls) if eltype(status.population) <: AbstractConstrainedSolution n = count(is_feasible, status.population) - @printf(io,"%12s %d / %d in final population\n", "feasibles:", n, length(status.population)) + @printf(io," %-16s %d / %d in final population\n", "Feasibles:", n, length(status.population)) end - @printf(io,"%12s %.4f s\n", "total time:", status.final_time - status.start_time) + @printf(io," %-16s %.4f s\n", "Total time:", status.overall_time) txt = status.stop ? termination_status_message(status) : "" - @printf(io,"%12s %s\n", "stop reason:", txt) - println(io, "+============================+") + @printf(io," %-16s %s\n", "Stop reason:", txt) end function Base.show(io::IO, ::MIME"text/html", status::State) println(io, "
")
-    show(io, "text/plain", status)
+    show(IOContext(io, :color => false), "text/plain", status)
     println(io, "
") end diff --git a/src/solutions/individual.jl b/src/solutions/individual.jl index ac84381c..3721ef99 100644 --- a/src/solutions/individual.jl +++ b/src/solutions/individual.jl @@ -9,6 +9,7 @@ mutable struct xf_solution{X} <: AbstractUnconstrainedSolution # Single Objectiv end const xf_indiv = xf_solution{Vector{Float64}} +include("utils.jl") # Single Objective Constrained mutable struct xfgh_solution{X} <: AbstractConstrainedSolution @@ -130,28 +131,28 @@ julia> population = [ Metaheuristics.create_child(rand(2), (randn(2), randn(2), ``` """ -function create_child(x, fResult::Float64;ε=0.0) - return xf_solution(x, fResult) +function create_child(x, fResult::T;ε=0.0) where T <: Real + return xf_solution(x, Float64(fResult)) end # constrained single objective function create_child(x, - fResult::Tuple{Float64,Array{Float64,1},Array{Float64,1}}; + fResult::Tuple{Real,Array{<:Real,1},Array{<:Real,1}}; ε=0.0 ) f, g, h = fResult isempty(g) && isempty(h) && (return xf_solution(x, f)) - return xfgh_indiv(x, f, g, h; ε=ε) + return xfgh_indiv(x, Float64(f), Float64.(g), Float64.(h); ε=ε) end # constrained multi-objective function create_child(x, - fResult::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}}; + fResult::Tuple{Vector{<:Real},Vector{<:Real},Vector{<:Real}}; ε = 0.0 ) f, g, h = fResult - return xFgh_indiv(x, f, g, h;ε=ε) + return xFgh_indiv(x, Float64.(f), Float64.(g), Float64.(h);ε=ε) end ##########################################################3 @@ -225,41 +226,35 @@ function create_child(X::AbstractMatrix, population end -##########################################################3 -function _gen_X(N, bounds::AbstractMatrix{T}) where T <: AbstractFloat - a = view(bounds, 1, :)' - b = view(bounds, 2, :)' - D = length(a) - a .+ (b - a) .* rand(eltype(bounds), N, D) +#= +function create_child(X::AbstractMatrix, fResult; ε=0.0) + error("""Objective function should return either a Vector or a Tuple depending on the problem. + Current output: $fResult of type $(typeof(fResult)) + Examples on parallelization at https://jmejia8.github.io/Metaheuristics.jl/stable/tutorials/parallelization""") end +=# -function _gen_X(N, bounds::AbstractMatrix{T}) where T <: Integer - a = view(bounds, 1, :)' - b = view(bounds, 2, :)' - D = length(a) - X = zeros(eltype(bounds), N, D) - for i in 1:D - X[:,i] = rand(a[i]:b[i], N) - end - return X +function create_child(x, fResult; ε = 0.0) + error(""" + Objective function should return either a numerical value or a Tuple depending on the problem. + Current output: `$fResult` of type `$(typeof(fResult))` + Examples at https://jmejia8.github.io/Metaheuristics.jl/stable/examples""") end -function _gen_X(N, bounds::AbstractMatrix{T}) where T <: Bool - rand(eltype(bounds), N, size(bounds, 2)) -end +##########################################################3 + +function generate_population(N::Int, problem, rng = default_rng_mh();ε=0.0, parallel_evaluation = false) -function generate_population(N::Int, problem;ε=0.0, parallel_evaluation = false) - X = _gen_X(N, problem.bounds) + space = problem.search_space if problem.parallel_evaluation + X = sample(RandomSampler(space; rng), N) return create_solutions(X, problem; ε=ε) end - population = [ create_solution(X[i,:], problem; ε=ε) for i in 1:N] - - return population + [create_solution(x, problem; ε=ε) for x in rand(rng, space, N)] end @@ -316,17 +311,19 @@ julia> population = [ Metaheuristics.create_child(rand(2), (randn(2), randn(2), generateChild(x, fx) = create_child(x, fx) function create_solution(x::AbstractVector, problem::Problem; ε=0.0) + xx = _fix_type(x, problem.search_space) problem.f_calls += 1 - return create_child(x, problem.f(x), ε=ε) + return create_child(xx, problem.f(xx), ε=ε) end function create_solutions(X::AbstractMatrix, problem::Problem; ε=0.0) - problem.f_calls += size(X,1) if problem.parallel_evaluation + problem.f_calls += size(X,1) + # FIXME: fix type of X return create_child(X, problem.f(X), ε=ε) end - [create_child(X[i,:], problem.f(X[i,:]), ε=ε) for i in 1:size(X,1)] + [create_solution(X[i,:], problem, ε=ε) for i in 1:size(X,1)] end # getters for the above structures diff --git a/src/solutions/utils.jl b/src/solutions/utils.jl new file mode 100644 index 00000000..54faf7d8 --- /dev/null +++ b/src/solutions/utils.jl @@ -0,0 +1,18 @@ +_fix_type(x::AbstractVector{T}, bounds::BoxConstrainedSpace{T}) where T = x + +function _fix_type( + x::AbstractVector{<:Real}, + bounds::BoxConstrainedSpace{T}, + ) where T <:Integer + round.(T, x) +end + +function _fix_type( + x::AbstractVector{<:Integer}, + bounds::BoxConstrainedSpace{T}, + ) where T <: AbstractFloat + T.(x) +end + +_fix_type(x, space) = x + diff --git a/src/common/stop.jl b/src/termination/budget.jl similarity index 92% rename from src/common/stop.jl rename to src/termination/budget.jl index 5fe34d29..a420b03e 100644 --- a/src/common/stop.jl +++ b/src/termination/budget.jl @@ -117,9 +117,3 @@ end -function default_stop_check(status, information, options) - call_limit_stop_check(status, information, options) || - iteration_stop_check(status, information, options) || - time_stop_check(status, information, options) || - accuracy_stop_check(status, information, options) -end diff --git a/src/termination/convergence.jl b/src/termination/convergence.jl new file mode 100644 index 00000000..402582c2 --- /dev/null +++ b/src/termination/convergence.jl @@ -0,0 +1,119 @@ +abstract type ConvergenceTermination <: AbstractTermination end + +include("moo_convergence.jl") + + + +function termination_status_message(criterion::ConvergenceTermination) + "Due to Convergence Termination criterion." +end + +function stop_check(status::State, criterion::ConvergenceTermination; report=true) + population = status.population + # nothing to do for empty populations + isempty(population) && return false + + feasible = is_feasible.(population) + + # check whether feasibility is ratio over 30% + count(feasible) / length(population) < 0.3 && return false + + stop = stop_check(population[feasible], criterion) + + + report && stop && (status.termination_status_code = criterion) + + stop +end + +""" + RelativeFunctionConvergence(;ftol) + +Check for a small relative difference between the function values of the +current population with the largest and smallest function values. +""" +Base.@kwdef struct RelativeFunctionConvergence <: ConvergenceTermination + ftol::Float64 = eps() +end + + +function stop_check(population::AbstractVector, criterion::RelativeFunctionConvergence) + + m, M = extrema(fvals(population)) + abs(M - m) / max(abs(M), eps()) <= criterion.ftol +end + +""" + SmallStandardDeviation(;ftol) + +Termination requires a small standard deviation of the function values of +the population. +""" +Base.@kwdef struct SmallStandardDeviation <: ConvergenceTermination + ftol::Float64 = 1e-6 +end + +function stop_check(population::AbstractVector, criterion::SmallStandardDeviation) + fx = fvals(population) + std(fx) <= criterion.ftol +end + +""" + AbsoluteFunctionConvergence(;ftol) + +Termination requires a small absolute difference between the function +values of the population with the largest and smallest function values. +""" +Base.@kwdef struct AbsoluteFunctionConvergence <: ConvergenceTermination + ftol::Float64 = 0.0 +end + +function stop_check(population::AbstractVector, criterion::AbsoluteFunctionConvergence) + m, M = extrema(fvals(population)) + abs(M - m) <= criterion.ftol +end + +""" + RelativeParameterConvergence(;xtol) + +Termination requires a small relative parameter difference between the +vertices with the largest and smallest function values. +""" +Base.@kwdef struct RelativeParameterConvergence <: ConvergenceTermination + xtol::Float64 = 1e-8 +end + +function stop_check(population::AbstractVector, criterion::RelativeParameterConvergence) + + fx = fvals(population) + xhi = population[argmax(fx)] |> get_position + xlo = population[argmin(fx)] |> get_position + # check for numerical representation + if eltype(xhi) <: Number + a = max(abs.(xlo), abs.(xhi)) + return maximum(abs.(xhi -xlo))/max(maximum(a), eps()) <= criterion.xtol + end + xhi == xlo +end + +struct CheckConvergence <: ConvergenceTermination + criteria::Vector +end + +function CheckConvergence(; + f_tol_rel = eps(), + f_tol_abs = 0.0, + x_tol = 1e-8, + criteria = [ + AbsoluteFunctionConvergence(f_tol_abs), + RelativeFunctionConvergence(f_tol_rel), + SmallStandardDeviation(), + RelativeParameterConvergence(x_tol), + ] + ) + CheckConvergence(criteria) +end + +function stop_check(population::AbstractVector, criteria::CheckConvergence) + all(stop_check(population, criterion) for criterion in criteria.criteria) +end diff --git a/src/termination/default.jl b/src/termination/default.jl new file mode 100644 index 00000000..56840254 --- /dev/null +++ b/src/termination/default.jl @@ -0,0 +1,20 @@ +function default_stop_check(status, information, options) + # budget termination + call_limit_stop_check(status, information, options) || + iteration_stop_check(status, information, options) || + time_stop_check(status, information, options) || + # accuracy + accuracy_stop_check(status, information, options) || + # convergence + stop_check(status, CheckConvergence( + f_tol_abs = options.f_tol, + f_tol_rel = options.f_tol_rel, + x_tol =options.x_tol)) #=|| + stop_check( + status, + CheckConvergence( + criteria=[RobustConvergence(ftol=options.f_tol)] + ) + ) + =# +end diff --git a/src/termination/moo_convergence.jl b/src/termination/moo_convergence.jl new file mode 100644 index 00000000..5b615e65 --- /dev/null +++ b/src/termination/moo_convergence.jl @@ -0,0 +1,51 @@ +abstract type MOOConvergenceTermination <: AbstractTermination end + + +Base.@kwdef mutable struct RobustConvergence <: ConvergenceTermination + ftol::Float64 = 1e-3 + period::Int = 10 + _last_population = [] + _iter::Int = 0 +end + +function _check_extremas(fa, fb, tol) + maximum(abs.(fb - fa)) <= tol +end + +function stop_check(f_current, f_last, criterion::RobustConvergence) + tol = criterion.ftol + # compute extrema + fmin = ideal(f_current) + fmax = nadir(f_current) + # normalization denominator + Δ = fmax - fmin + Δ[Δ .≈ 0] .= one(eltype(Δ)) + + # if extrema are not in tolerance, then don't stop + !_check_extremas(fmin ./ Δ, ideal(f_last) ./ Δ, tol) && return false + !_check_extremas(fmax ./ Δ, nadir(f_last) ./ Δ, tol) && return false + + + + # normalize fronts + f_current = (f_current .- fmin') ./ Δ' + f_last = (f_last .- fmin') ./ Δ' + # stop if extrema and now IGD is under tolerance + indicator = PerformanceIndicators.igd(f_last, f_current) + indicator <= tol +end + +function stop_check(population::AbstractVector, criterion::RobustConvergence) + last_population = criterion._last_population + if isempty(last_population) + criterion._last_population = deepcopy(population) + end + + criterion._iter += 1 + criterion._iter % criterion.period != 0 && return false + + stop = stop_check(fvals(population), fvals(last_population), criterion) + criterion._last_population = deepcopy(population) + stop +end + diff --git a/src/termination/stop_status_codes.jl b/src/termination/stop_status_codes.jl new file mode 100644 index 00000000..eac247e5 --- /dev/null +++ b/src/termination/stop_status_codes.jl @@ -0,0 +1,57 @@ +struct IterationLimit <: AbstractTermination end +struct TimeLimit <: AbstractTermination end +struct EvaluationsLimit <: AbstractTermination end +struct AccuracyLimit <: AbstractTermination end +struct OtherLimit <: AbstractTermination end +struct ObjectiveVarianceLimit <: AbstractTermination end +struct ObjectiveDifferenceLimit <: AbstractTermination end +struct UnknownStopReason <: AbstractTermination end + + +const ITERATION_LIMIT = IterationLimit() +const TIME_LIMIT = TimeLimit() +const EVALUATIONS_LIMIT = EvaluationsLimit() +const ACCURACY_LIMIT = AccuracyLimit() +const OTHER_LIMIT = OtherLimit() +const UNKNOWN_STOP_REASON = UnknownStopReason() +const OBJECTIVE_VARIANCE_LIMIT = ObjectiveVarianceLimit() +const OBJECTIVE_DIFFERENCE_LIMIT= ObjectiveDifferenceLimit() + +""" + termination_status_message(status) + +Return a string of the message related to the `status`. + +See [`TerminationStatusCode`](@ref) + +### Example: + +```julia-repl +julia> termination_status_message(Metaheuristics.ITERATION_LIMIT) +"Maximum number of iterations exceeded." + +julia> termination_status_message(optimize(f, bounds)) +"Maximum number of iterations exceeded." + +julia> termination_status_message(ECA()) +"Unknown stop reason." +``` +""" +function termination_status_message(status_code::TerminationStatusCode) + codes = Dict( + ITERATION_LIMIT => "Maximum number of iterations exceeded.", + TIME_LIMIT => "Maximum time exceeded.", + EVALUATIONS_LIMIT => "Maximum objective function calls exceeded.", + ACCURACY_LIMIT => "The desired accuracy was obtained.", + OBJECTIVE_VARIANCE_LIMIT => "Small variance of the objective function.", + OBJECTIVE_DIFFERENCE_LIMIT =>"Small difference of objective function values.", + OTHER_LIMIT => "Other stopping criteria.", + UNKNOWN_STOP_REASON => "Unknown stop reason." + ) + try + return codes[status_code] + catch + throw("Illegal value $status_code") + end +end + diff --git a/src/termination/termination.jl b/src/termination/termination.jl new file mode 100644 index 00000000..9e699ca5 --- /dev/null +++ b/src/termination/termination.jl @@ -0,0 +1,27 @@ +include("stop_status_codes.jl") +include("convergence.jl") +include("budget.jl") +include("default.jl") + +""" + Termination(checkall, checkany) + +Define a termination criterion. This criterion will suggest a stop condition +if all conditions in `checkall` are satisfied or at least one in `checkany` +is fulfilled. +""" +Base.@kwdef struct Termination <: AbstractTermination + checkall::Vector = Any[] + checkany::Vector = Any[] +end + +function stop_check(status, termination::Termination) + checkall = termination.checkall + # all criteria are satisfied + !isempty(checkall) && return all(stop_check(status, c) for c in checkall) + + # check if any other does + any(stop_check(status, c) for c in termination.checkany) +end + + diff --git a/test/multi-objective.jl b/test/multi-objective.jl index c623ef38..513c542c 100644 --- a/test/multi-objective.jl +++ b/test/multi-objective.jl @@ -47,7 +47,7 @@ end @test Metaheuristics.PerformanceIndicators.igd(result.population, pf) <= 0.5 @test Metaheuristics.PerformanceIndicators.spacing(result) < 0.5 @test Metaheuristics.PerformanceIndicators.covering(pf, result.population) <= 1.0 - @test Metaheuristics.PerformanceIndicators.covering(result, result) ≈ 0.0 + @test Metaheuristics.PerformanceIndicators.covering(result, result, verbose=false) ≈ 0.0 # number of function evaluations should be reported correctly @test f_calls == result.f_calls