diff --git a/scripts/ga_testing.jl b/scripts/ga_testing.jl new file mode 100644 index 0000000..c59cf73 --- /dev/null +++ b/scripts/ga_testing.jl @@ -0,0 +1,142 @@ + +using GlobalOptimization +using BenchmarkTools +using Random +using Distributions +#using BlackBoxOptim +#using LoopVectorization +#using PaddedViews +#using StaticArrays +using Profile +#using JET +using Infiltrator +#Random.seed!(1234) + +# Schwefel Function +function schaffer(x) + obj = 0.5 + (sin(x[1]^2 + x[2]^2)^2 - 0.5) / (1 + 0.001 * (x[1]^2 + x[2]^2))^2 + return obj, 0.0 +end + +function waveDrop(x) + obj = -(1 + cos(12 * sqrt(x[1]^2 + x[2]^2))) / (0.5 * (x[1]^2 + x[2]^2) + 2.0) + return obj, 0.0 +end + +@inline function layeb_1(x) + obj = 0.0 + @fastmath for val in x + xm1sq = (val - 1)^2 + obj += 10000.0 * sqrt(abs(exp(xm1sq) - 1.0)) + end + return obj, 0.0 +end + +function rastrigin(x; A=10) + obj = A * length(x) + for val in x + obj += val^2 - A * cos(2 * pi * val) + end + #sleep(rand()*2e-6) + return obj +end + +# Setup Problem +N = 7 # numdims +ss = ContinuousRectangularSearchSpace([-100.0 for i in 1:N], [100.0 for i in 1:N]) +prob = GlobalOptimization.OptimizationProblem(rastrigin, ss) + +# Instantiate DE +mutation_probability = 0.05; +mutation_delta = 2 +mutation_strategy = RandomElementwiseMutation( + mutation_delta; dist=Normal(mutation_probability, 0.001) +) +crossover_strategy = BLXAlphaCrossover(.8, 0.5) +elitism_strategy = NoElitism() +selection_strategy = TournamentSelection() +ga = GA( + prob; + eval_method=SerialFunctionEvaluation(), + num_candidates=100, + max_iterations=1000, + max_stall_iterations=100, + mutation_params=mutation_strategy, + selection_params=selection_strategy, + crossover_params=crossover_strategy, + elitism_params=elitism_strategy, + show_trace=Val(false), +) + +res = @btime optimize!(ga) + +# run 100 times and get avg +println("beginning batch test") +fbests = Float64[] +for i = 1:100 + #println("i=$i") + gai = GA( + prob; + eval_method=SerialFunctionEvaluation(), + num_candidates=100, + max_iterations=1000, + max_stall_iterations=100, + mutation_params=mutation_strategy, + selection_params=selection_strategy, + crossover_params=crossover_strategy, + elitism_params=elitism_strategy, + show_trace=Val(false), + ) + + resi = optimize!(gai) + push!(fbests, resi.fbest) +end +println("Avg Value: $(sum(fbests)/length(fbests))") +#iters_per_solve = map(i->optimize!(deepcopy(de)).iters, 1:100); + +# bb_res = bboptimize( +# rastrigin; +# Method = :adaptive_de_rand_1_bin_radiuslimited, +# PopulationSize = 100, +# SearchRange = (-5.0, 5.0), +# NumDimensions = N, +# TraceMode = :compact, +# TraceInterval = 0.001, +# #MaxSteps = 1000, +# ) + +# Instantiate PSO +# spso = SerialPSO(prob; max_time = 20.0) +#tpso = ThreadedPSO(prob; max_time = 20.0) +#ppso = PolyesterPSO(prob; max_time = 20.0) + +# #res = optimize!(spso) +# res = optimize!(spso); display(res) +# res = optimize!(tpso); display(res) +# res = optimize!(ppso); display(res) + +# ======== BENCHMARKING +# sres = @benchmark optimize!(_pso) setup=(_pso = SerialPSO(prob; max_iterations=20)) +# tres = @benchmark optimize!(_pso) setup=(_pso = ThreadedPSO(prob; max_iterations=20)) +# pres = @benchmark optimize!(_pso) setup=(_pso = PolyesterPSO(prob; max_iterations=20)) +# display(sres) +# display(tres) +# display(pres) +# GlobalOptimization.initialize!(spso) +# GlobalOptimization.update_velocity!(spso.swarm, spso.cache, 10, 0.5, 0.49, 0.49) +# GlobalOptimization.step!(spso.swarm) +# GlobalOptimization.enforce_bounds!(spso.swarm, ss) + +# go = @benchmark GlobalOptimization.evaluate_fitness!($spso.swarm, $spso.evaluator) +# display(go) + +# ======== ALLOCATION TRACKING +# pso1 = SerialPSO(prob) +# pso2 = SerialPSO(prob) +# optimize!(pso1) +# Profile.clear_malloc_data() +# optimize!(pso2) + +# ======== TYPES +#@report_call GlobalOptimization.optimize!(spso) +#report_package(GlobalOptimization) diff --git a/scripts/pso_testing.jl b/scripts/pso_testing.jl index 01c17f6..974fb4b 100644 --- a/scripts/pso_testing.jl +++ b/scripts/pso_testing.jl @@ -43,6 +43,11 @@ N = 10 ss = ContinuousRectangularSearchSpace([-5.12 for i in 1:N], [5.12 for i in 1:N]) prob = OptimizationProblem(layeb_1, ss) +# Rast +N = 7 +ss = ContinuousRectangularSearchSpace([-100 for i in 1:N], [100 for i in 1:N]) +prob = OptimizationProblem(rastrigin, ss) + # Instantiate PSO spso = PSO( prob; @@ -54,7 +59,7 @@ spso = PSO( tpso = PSO(prob; eval_method=ThreadedFunctionEvaluation(), max_time=20.0) ppso = PSO(prob; eval_method=PolyesterFunctionEvaluation(), max_time=20.0) -Random.seed!(1234) +#Random.seed!(1234) res = optimize!(spso) # res = optimize!(spso); display(res) # res = optimize!(tpso); display(res) diff --git a/src/GA/GA.jl b/src/GA/GA.jl new file mode 100644 index 0000000..a4cca92 --- /dev/null +++ b/src/GA/GA.jl @@ -0,0 +1,244 @@ +""" +GAOptions + +Options for the Genetic Algorithm (GA) algorithms. + +# Fields: +- `general<:GeneralOptions`: The general options. +- `pop_init_method<:AbstractPopulationInitialization`: The population initialization method. +- `selection_params<:AbstractGASelectionParameters`: The selection parameters +- `crossover_params<:AbstractGACrossoverParameters`: The crossover strategy parameters. +- `mutation_params<:AbstractGAMutationParameters`: The mutation strategy parameters. +- `initial_space<:Union{Nothing,ContinuousRectangularSearchSpace}`: The initial space to initialize the population. +""" +struct GAOptions{ + SP<:AbstractGASelectionParameters, + CP<:AbstractGACrossoverParameters, + MP<:AbstractGAMutationParameters, + EP<:AbstractGAElitismParameters, + ISS<:Union{Nothing,ContinuousRectangularSearchSpace}, + PI<:AbstractPopulationInitialization, + GO<:GeneralOptions, +} <: AbstractAlgorithmSpecificOptions + # The general options + general::GO + + # The Population initialization method + pop_init_method::PI + + # Selection params + selection_params::SP + + # Crossover params + crossover_params::CP + + # Mutation strategy parameters + mutation_params::MP + + elitism_params::EP + + # Initial space to initialize the population + initial_space::ISS + + """ + GAOptions(args...) + + Construct the Genetic Algorithm (GA) algorithms options. + + # Arguments + - `general<:GeneralOptions`: The general options. + - `pim<:AbstractPopulationInitialization`: The population initialization method. + - `selection<:AbstractGASelectionParameters`: The selection parameters + - `crossover<:AbstractGACrossoverParameters`: The crossover strategy parameters. + - `mutation<:AbstractGAMutationParameters`: The mutation strategy parameters. + - `elitism<:AbstractGAElitismParameters`: Elitism Parameters + - `initial_space<:Union{Nothing,ContinuousRectangularSearchSpace}`: The initial space to initialize the population. + """ + function GAOptions( + general::GO, pim::PI, selection::SP, crossover::CP, mutation::MP, elitism::EP, initial_space::ISS + ) where {SP<:AbstractGASelectionParameters, MP<:AbstractGAMutationParameters,CP<:AbstractGACrossoverParameters,EP<:AbstractGAElitismParameters,GO,PI,ISS} + return new{SP,CP,MP,EP,ISS,PI,GO}(general, pim, selection, crossover, mutation, elitism, initial_space) + end +end + + +""" + GA + +The genetic algorithm +""" +struct GA{ + SP<:AbstractGASelectionParameters, + CP<:AbstractGACrossoverParameters, + MP<:AbstractGAMutationParameters, + EP<:AbstractGAElitismParameters, + T<:AbstractFloat, + E<:BatchEvaluator, + IBSS, + PI<:AbstractPopulationInitialization, + GO<:GeneralOptions, +} <: AbstractPopulationBasedOptimizer + # The GA algorithm options + options::GAOptions{SP,CP,MP,EP,IBSS,PI,GO} + + # The GA evaluator + evaluator::E + + # The population + population::GAPopulation{T} + + # The GA cache + cache::MinimalPopulationBasedOptimizerCache{T} +end + + + +""" + GA(prob::AbstractProblem{has_penalty,SS}; kwargs...) + +Construct a serial Genetic Algorithm (GA) algorithm with the given options. + +# Arguments +- `prob::AbstractProblem{has_penalty,SS}`: The problem to solve. + +# Keyword Arguments +- `eval_method::AbstractFunctionEvaluationMethod=SerialFunctionEvaluation()`: The method to use for evaluating the objective function. +- `num_candidates::Integer=100`: The number of candidates in the population. +- `population_initialization::AbstractPopulationInitialization=UniformInitialization()`: The population initialization method. +- `selection_params::SP=StochasticUniversalSampling(100)`: The selection parameters. +- `crossover_params::CP=BLXAlphaCrossover(0.8, 0.5)`: The crossover strategy parameters. +- `mutation_params::MP=RandomElementwiseMutation(0.05, 0.01)`: The mutation strategy parameters. +- `elitism_params::EP=SimpleElitism(0.03)`: the elitism parameters +- `initial_space::Union{Nothing,ContinuousRectangularSearchSpace}=nothing`: The initial bounds for the search space. +- `max_iterations::Integer=1000`: The maximum number of iterations. +- `function_tolerance::Real=1e-6`: The function tolerance (stall-based stopping criteria). +- `max_stall_time::Real=60.0`: The maximum stall time (in seconds). +- `max_stall_iterations::Integer=100`: The maximum number of stall iterations. +- `max_time::Real=60.0`: The maximum time (in seconds) to allow for optimization. +- `min_cost::Real=(-Inf)`: The minimum cost to allow for optimization. +- `function_value_check::Union{Val{false},Val{true}}=Val(true)`: Whether to check the function value + for bad values (i.e., Inf or NaN). +- `show_trace::Union{Val{false},Val{true}}=Val(false)`: Whether to show the trace. +- `save_trace::Union{Val{false},Val{true}}=Val(false)`: Whether to save the trace. +- `save_file::String="trace.txt"`: The file to save the trace to. +- `trace_level::TraceLevel=TraceMinimal(1)`: The trace level to use. +""" +function GA( + prob::AbstractProblem{has_penalty,SS}; + eval_method::AbstractFunctionEvaluationMethod=SerialFunctionEvaluation(), + num_candidates::Integer=100, + population_initialization::AbstractPopulationInitialization=UniformInitialization(), + selection_params::SP=StochasticUniversalSampling(100), + crossover_params::CP=BLXAlphaCrossover(0.8, 0.5), + mutation_params::MP=RandomElementwiseMutation(0.05, 0.01), + elitism_params::EP=NoElitism(), + initial_space::Union{Nothing,ContinuousRectangularSearchSpace}=nothing, + max_iterations::Integer=1000, + function_tolerance::Real=1e-6, + max_stall_time::Real=60.0, + max_stall_iterations::Integer=100, + max_time::Real=60.0, + min_cost::Real=(-Inf), + function_value_check::Union{Val{false},Val{true}}=Val(true), + show_trace::Union{Val{false},Val{true}}=Val(false), + save_trace::Union{Val{false},Val{true}}=Val(false), + save_file::String="trace.txt", + trace_level::TraceLevel=TraceMinimal(1), +) where { + SP<:AbstractGASelectionParameters, + MP<:AbstractGAMutationParameters, + CP<:AbstractGACrossoverParameters, + EP<:AbstractGAElitismParameters, + T<:AbstractFloat, + SS<:ContinuousRectangularSearchSpace{T}, + has_penalty, +} + + # Construct options + options = GAOptions( + GeneralOptions( + GlobalOptimizationTrace(show_trace, save_trace, save_file, trace_level), + function_value_check, + min_cost, + max_time, + max_iterations, + function_tolerance, + max_stall_time, + max_stall_iterations, + ), + population_initialization, + selection_params, + crossover_params, + mutation_params, + elitism_params, + intersection(search_space(prob), initial_space), + ) + + # Construct evaluator + return GA( + options, + construct_batch_evaluator(eval_method, prob), + GAPopulation{T}(num_candidates, num_dims(prob)), + MinimalPopulationBasedOptimizerCache{T}(num_dims(prob)), + ) +end + +# ===== AbstractPopBasedOptimizer interface +get_population(opt::GA) = opt.population.current_generation + +# Initialize parameters and pop and evaluate pop once +function initialize!(opt::GA) + # Unpack DE + @unpack options, evaluator, population, cache = opt + @unpack pop_init_method, mutation_params, crossover_params, selection_params = options + + # Initialize population + initialize!(population, pop_init_method, options.initial_space) + + # Initialize selection, crossover, and mutation + initialize!(mutation_params) + initialize!(crossover_params) + initialize!(selection_params, population) + + + # Handle fitness + initialize_fitness!(population, evaluator) + check_fitness!(population.current_generation, get_function_value_check(options)) + + # Initialize cache + update_global_best!(opt) + initialize!(cache) + + return nothing +end + +function step!(opt::GA) + # Unpack GA + @unpack options, evaluator, population, cache = opt + @unpack mutation_params, crossover_params, selection_params, elitism_params = options + search_space = evaluator.prob.ss + + # Perform selection + selection!(population, selection_params) + + # Perform crossover + crossover!(population, crossover_params) + + # Perform mutation + mutate!(population, mutation_params, search_space) + + # Evaluate fitness + evaluate_children!(population, evaluator) + check_fitness!(population.children, get_function_value_check(options)) + + set_next_generation!(population, elitism_params) + + # Update global best + improved = update_global_best!(opt) + + # Adapt mutation and crossover parameters if necessary + adapt!(mutation_params, improved) + adapt!(crossover_params, improved) + + return nothing +end \ No newline at end of file diff --git a/src/GA/Population.jl b/src/GA/Population.jl new file mode 100644 index 0000000..f9b8d31 --- /dev/null +++ b/src/GA/Population.jl @@ -0,0 +1,91 @@ +""" + GABasePopulation{T<:AbstractFloat} <: AbstractPopulation + +Base representation of population for Genetic Algorithms. This allcoates storage for a single set of candidates. +Multiple of these are used within the GA to represent the population, mating pool, and children. +""" +struct GABasePopulation{T<:AbstractFloat} <: AbstractPopulation{T} + + candidates::Vector{Vector{T}} + candidates_fitness::Vector{T} + + function GABasePopulation{T}(num_candidates::Integer, num_dims::Integer) where {T} + return new{T}( + [zeros(T, num_dims) for _ in 1:num_candidates], zeros(T, num_candidates) + ) + end + +end + +""" + GAPopulation{T<:AbstractFloat} +Full representation for GA population, with space for candidates and mating pool +""" +struct GAPopulation{T<:AbstractFloat} + # The current population of candidates + current_generation::GABasePopulation{T} + + # The mating pool, to be populated in the selection step + mating_pool::GABasePopulation{T} + + # children: generated and operated on in the crossover step + children::GABasePopulation{T} + + # Vectors of integers to store indexes for elitism/selection purposes + parent_idx_array::Vector{Int} # parents + child_idx_array::Vector{Int} # children + + function GAPopulation{T}(num_candidates::Integer, num_dims::Integer) where {T} + num_dims > 0 || throw(ArgumentError("num_dims must be greater than 0.")) + num_candidates > 0 || throw(ArgumentError("num_candidates must be greater than 0.")) + return new{T}( + GABasePopulation{T}(num_candidates, num_dims), + GABasePopulation{T}(num_candidates, num_dims), + GABasePopulation{T}(num_candidates, num_dims), + Vector{Int}(zeros(num_candidates)), + Vector{Int}(zeros(num_candidates)) + ) + end +end + +""" + function GAPopulation(num_candidates::Integer, num_dims::Integer) + Creates a `GAPopulation` with `num_candidates`, each having `num_dims`. +""" +function GAPopulation(num_candidates::Integer, num_dims::Integer) + return GAPopulation{Float64}(num_candidates, num_dims) +end + +Base.length(population::GAPopulation) = length(population.current_generation) + +function initialize!( + population::GAPopulation{T}, + pop_init_method::AbstractPopulationInitialization, + search_space::ContinuousRectangularSearchSpace{T}, +) where {T} + # Unpack population + @unpack candidates = population.current_generation + @unpack dim_min, dim_max = search_space + + initialize_population_vector!(candidates, dim_min, dim_max, pop_init_method) + + return nothing +end + +function initialize_fitness!( + population::GAPopulation{T}, evaluator::BatchEvaluator +) where {T} + # Evaluate the cost function for each candidate + evaluate!(population.current_generation, evaluator) + return nothing +end + +function evaluate_fitness!(population, evaluator) + evaluate!(population.current_generation, evaluator) + return nothing +end + +function evaluate_children!(population, evaluator) + evaluate!(population.children, evaluator) + return nothing +end \ No newline at end of file diff --git a/src/GA/crossover.jl b/src/GA/crossover.jl new file mode 100644 index 0000000..7a5f64d --- /dev/null +++ b/src/GA/crossover.jl @@ -0,0 +1,134 @@ + +""" + AbstractGACrossoverParameters{AS<:AbstractAdaptationStrategy} + +Abstract type for the crossover parameters used within the Genetic Algorithm (GA). + +Subtypes of this abstract type should define the following methods: +- `get_parameter(params::AbstractGACrossoverParameters, i)`: Returns the crossover parameter for the `i`-th candidate. +- `initialize!(params::AbstractGACrossoverParameters, num_dims, population_size)`: Initializes the crossover parameters. +- `adapt!(params::AbstractGACrossoverParameters, improved, global_best_improved)`: Adapts the crossover parameters based on the improvement status of the candidates. +- `crossover!(population::GAPopulation, crossover_params, search_space)`: Performs the crossover operation on the mating pool using the specified crossover parameters. +""" + +abstract type AbstractGACrossoverParameters{AS<:AbstractAdaptationStrategy} end + +function adapt!( + params::AbstractGACrossoverParameters{NoAdaptation}, global_best_improved +) + return nothing +end + +function adapt!( + params::AbstractGACrossoverParameters{RandomAdaptation}, global_best_improved +) + + if !global_best_improved + params.CR = one_clamped_rand(params.dist) + end + + return nothing + +end + +""" + BLXAlphaCrossover{AS<:AbstractAdaptationStrategy} + +The BLX-α crossover mating strategy for GA. +""" +mutable struct BLXAlphaCrossover{AS<:AbstractAdaptationStrategy, D} <: AbstractGACrossoverParameters{AS} + + CR::Float64 + dist::D # in case we desire an adaptive parameter + alpha::Float64 # BKX-α parameter + + + """ + function BLXAlphaCrossover(CR::Float64) + Creates non-adaptive BLX-α crossover parameters with a static crossover ratio CR applied to the whole population. + """ + function BLXAlphaCrossover(CR::Float64, alpha) + return new{NoAdaptation,Nothing}(CR, nothing, alpha) + end + + + """ + function BLXAlphaCrossover(; dist=default_ga_crossover_dist) + Creates an adaptive BLX-α crossover parameters where `CR` is determined by `dist`. + """ + function BLXAlphaCrossover(alpha; dist=default_ga_crossover_dist) + return new{RandomAdaptation, typeof(dist)}(0.0, dist, alpha) # intialize CR as 0.0, will be changed in initialize! + end + +end + +function initialize!(params::BLXAlphaCrossover{NoAdaptation}) # For no adaptation, no initialization is needed + return nothing +end + +function initialize!(params::BLXAlphaCrossover{RandomAdaptation}) + params.CR = rand(params.dist) + return nothing +end + + +""" + function crossover!(population::GAPopulation, crossover_params::BLXAlphaCrossover, search space) + + Performs the BLX-α crossover on the `population` +""" +function crossover!( + population::GAPopulation, + crossover_params::BLXAlphaCrossover +) + + # Unpack the 3 sub populations: current generation, mating pool, children. We will be populating the children in this function + @unpack current_generation, mating_pool, children = population + + # Point to the candidates of each + current_generation = current_generation.candidates + mating_pool = mating_pool.candidates + children = children.candidates + + alpha = crossover_params.alpha + m = length(current_generation[1]) + + # first, need to see how many pairs we can make + pop_len = length(population) + + extra_candidates = mod(pop_len, 2) + pop_is_odd = extra_candidates == 1 # if pop is odd (why) we'll have one sad, lonely candidate we'll need to handle + + # shuffle the MP + shuffle!(mating_pool) + + # shuffle was random, so we can iterate over pairs + @inbounds for i=1:2:pop_len-1-extra_candidates # if any extra candidates, we'll leave it off at the end and add it to children + + if rand()current_generation.candidates_fitness[i]) + + # Assign best performer to mating pool + # OPTIONAL: later, include some probability based selection here + mating_pool.candidates[t] = current_generation.candidates[participants[1]] + mating_pool.candidates_fitness[t] = current_generation.candidates_fitness[participants[1]] + end + + return nothing + +end \ No newline at end of file diff --git a/src/GA/util.jl b/src/GA/util.jl new file mode 100644 index 0000000..10db37f --- /dev/null +++ b/src/GA/util.jl @@ -0,0 +1,4 @@ +# a basic default truncated distribution. not advised by anything except centered around 0.8, which seems to be best based on +# the MAE552 notes. The truncated points are where the distribution approaches zero, so this should be fine. +const default_ga_crossover_dist = Normal(0.9, 0.1) +const default_ga_mutation_dist = Normal(0.025, 0.005) \ No newline at end of file diff --git a/src/GlobalOptimization.jl b/src/GlobalOptimization.jl index 555f0da..904b4f8 100644 --- a/src/GlobalOptimization.jl +++ b/src/GlobalOptimization.jl @@ -1,7 +1,7 @@ module GlobalOptimization using ChunkSplitters: chunks, ChunkSplitters, RoundRobin -using Distributions: Cauchy, Laplace, MixtureModel +using Distributions: Cauchy, Laplace, MixtureModel, Normal, sample using LatinHypercubeSampling: scaleLHC, LHCoptim using LinearAlgebra: dot, eigen!, mul!, tril! using Polyester: @batch @@ -48,6 +48,15 @@ include("MBH/Hopper.jl") include("MBH/LocalSearch.jl") include("MBH/MBH.jl") +# GA +include("GA/util.jl") +include("GA/Population.jl") +include("GA/mutation.jl") +include("GA/crossover.jl") +include("GA/selection.jl") +include("GA/elitism.jl") +include("GA/GA.jl") + export TraceMinimal, TraceDetailed, TraceAll export ContinuousRectangularSearchSpace @@ -56,10 +65,12 @@ export LatinHypercubeInitialization export OptimizationProblem, NonlinearProblem, NonlinearLeastSquaresProblem export optimize! -export PSO, DE, MBH +export PSO, DE, MBH, GA +# PSO exports export MATLABVelocityUpdate, CSRNVelocityUpdate +# DE exports export SimpleSelector, RadiusLimitedSelector, RandomSubsetSelector export Rand1, Rand2, Best1, Best2, CurrentToBest1, CurrentToBest2 export CurrentToRand1, CurrentToRand2, RandToBest1, RandToBest2, Unified @@ -67,10 +78,15 @@ export MutationParameters, SelfMutationParameters export SelfBinomialCrossoverParameters, BinomialCrossoverParameters export CovarianceTransformation, UncorrelatedCovarianceTransformation +# MBH exports export SingleHopper, MCH export MBHStaticDistribution, MBHAdaptiveDistribution export LocalStochasticSearch, UserLocalSearch +# GA exports +export StochasticUniversalSampling, TournamentSelection, BLXAlphaCrossover, RandomElementwiseMutation +export NoElitism, SimpleElitism + # Handle extension symbols we want to export # NOTE: I really don't like this solution, but it seems to be the best option for now... abstract type BOBYQALocalSearch{T} <: DerivativeBasedLocalSearch{T} end