allow better control of RNG (part 2)

- use `rng` option in DE & ES
- updates tests to use `rng` option

Author: wildart

README.md

@@ -26,7 +26,7 @@ pkg> add Evolutionary
 - Genetic Programming (TreeGP)
 
 ## Operators
-a
+
 - Mutations
   - ES
     - (an)isotropic Gaussian

src/de.jl

@@ -48,6 +48,7 @@ function update_state!(objfun, constraints, state, population::AbstractVector{IT
     Np = method.populationSize
     n = method.n
     F = method.F
+    rng = options.rng
 
     offspring = Array{IT}(undef, Np)
 
@@ -61,12 +62,12 @@ function update_state!(objfun, constraints, state, population::AbstractVector{IT
         offspring[i] = copy(base)
         # println("$i => base:", offspring[i])
 
-        targets = randexcl(1:Np, [i], 2*n)
+        targets = randexcl(rng, 1:Np, [i], 2*n)
         offspring[i] = differentiation(offspring[i], @view population[targets]; F=F)
         # println("$i => mutated:", offspring[i], ", targets:", targets)
 
         # recombination
-        offspring[i], _ = method.recombination(offspring[i], base)
+        offspring[i], _ = method.recombination(offspring[i], base, rng=rng)
         # println("$i => recombined:", offspring[i])
     end
 

src/es.jl

@@ -6,8 +6,8 @@ The constructor takes following keyword arguments:
 - `initStrategy`: an initial strategy description, (default: empty)
 - `recombination`: ES recombination function for population (default: `first`), see [Crossover](@ref)
 - `srecombination`: ES recombination function for strategies (default: `first`), see [Crossover](@ref)
-- `mutation`: [Mutation](@ref) function for population (default: `first`)
-- `smutation`: [Mutation](@ref) function for strategies (default: `identity`)
+- `mutation`: [Mutation](@ref) function for population (default: [`nop`](@ref))
+- `smutation`: [Mutation](@ref) function for strategies (default: [`nop`](@ref))
 - `μ`/`mu`: the number of parents
 - `ρ`/`rho`: the mixing number, ρ ≤ μ, (i.e., the number of parents involved in the procreation of an offspring)
 - `λ`/`lambda`: the number of offspring
@@ -27,8 +27,8 @@ struct ES <: AbstractOptimizer
     ES(; initStrategy::AbstractStrategy = NoStrategy(),
         recombination::Function = first,
         srecombination::Function = first,
-        mutation::Function = ((r,s) -> r),
-        smutation::Function = identity,
+        mutation::Function = nop,
+        smutation::Function = nop,
         μ::Integer = 1,
         mu::Integer = μ,
         ρ::Integer = μ,
@@ -74,6 +74,7 @@ end
 function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::ES, options, itr) where {IT}
     @unpack initStrategy,recombination,srecombination,mutation,smutation,μ,ρ,λ,selection = method
     evaltype = options.parallelization
+    rng = options.rng
 
     @assert ρ <= μ "Number of parents involved in the procreation of an offspring should be no more then total number of parents"
     if selection == :comma
@@ -86,21 +87,21 @@ function update_state!(objfun, constraints, state, population::AbstractVector{IT
     for i in 1:λ
         # Recombine the ρ selected parents to form a recombinant individual
         if ρ == 1
-            j = rand(1:μ)
+            j = rand(rng, 1:μ)
             recombinantStrategy = state.strategies[j]
             recombinant = copy(population[j])
         else
-            idx = randperm(μ)[1:ρ]
+            idx = randperm(rng, μ)[1:ρ]
             recombinantStrategy = srecombination(state.strategies[idx])
-            recombinant = recombination(population[idx])
+            recombinant = recombination(population[idx]; rng=rng)
         end
 
         # Mutate the strategy parameter set of the recombinant
-        stgoff[i] = smutation(recombinantStrategy)
+        stgoff[i] = smutation(recombinantStrategy; rng=rng)
 
         # Mutate the objective parameter set of the recombinant using the mutated strategy parameter set
         # to control the statistical properties of the object parameter mutation
-        off = mutation(recombinant, stgoff[i])
+        off = mutation(recombinant, stgoff[i]; rng=rng)
 
         # Apply constraints
         offspring[i] = apply!(constraints, off)
@@ -136,3 +137,4 @@ function update_state!(objfun, constraints, state, population::AbstractVector{IT
 
     return false
 end
+

src/ga.jl

@@ -58,54 +58,72 @@ function initial_state(method::GA, options, objfun, population)
     return GAState(N, eliteSize, minfit, fitness, copy(population[fitidx]))
 end
 
-function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::GA, options, itr) where {IT}
-    @unpack populationSize,crossoverRate,mutationRate,ɛ,selection,crossover,mutation = method
+function update_state!(objfun, constraints, state, parents::AbstractVector{IT}, method::GA, options, itr) where {IT}
+    populationSize = method.populationSize
     evaltype = options.parallelization
     rng = options.rng
-    offspring = similar(population)
+    offspring = similar(parents)
 
-    # Select offspring
-    selected = selection(state.fitpop, populationSize, rng=rng)
+    # select offspring
+    selected = method.selection(state.fitpop, populationSize, rng=rng)
 
-    # Perform mating
-    offidx = randperm(rng, populationSize)
+    # perform mating
     offspringSize = populationSize - state.eliteSize
-    for i in 1:2:offspringSize
-        j = (i == offspringSize) ? i-1 : i+1
-        if rand(rng) < crossoverRate
-            offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]], rng=rng)
-        else
-            offspring[i], offspring[j] = population[selected[i]], population[selected[j]]
-        end
-    end
+    recombine!(offspring, parents, selected, method, offspringSize)
 
     # Elitism (copy population individuals before they pass to the offspring & get mutated)
     fitidxs = sortperm(state.fitpop)
     for i in 1:state.eliteSize
         subs = offspringSize+i
-        offspring[subs] = copy(population[fitidxs[i]])
+        offspring[subs] = copy(parents[fitidxs[i]])
     end
 
-    # Perform mutation
-    for i in 1:offspringSize
-        if rand(rng) < mutationRate
-            mutation(offspring[i], rng=rng)
-        end
-    end
+    # perform mutation
+    mutate!(offspring, method, constraints, rng=rng)
 
-    # Create new generation & evaluate it
-    for i in 1:populationSize
-        population[i] = apply!(constraints, offspring[i])
-    end
     # calculate fitness of the population
-    value!(objfun, state.fitpop, population)
-    # apply penalty to fitness
-    penalty!(state.fitpop, constraints, population)
+    evaluate!(objfun, state.fitpop, offspring, constraints)
 
-    # find the best individual
+    # select the best individual
     minfit, fitidx = findmin(state.fitpop)
-    state.fittest = population[fitidx]
+    state.fittest = parents[fitidx]
     state.fitness = state.fitpop[fitidx]
+    
+    # replace population
+    parents .= offspring
 
     return false
 end
+
+function recombine!(offspring, parents, selected, method, n=length(selected);
+                    rng::AbstractRNG=Random.GLOBAL_RNG)
+    mates = ((i,i == n ? i-1 : i+1) for i in 1:2:n)
+    for (i,j) in mates
+        p1, p2 = parents[selected[i]], parents[selected[j]]
+        if rand(rng) < method.crossoverRate
+            offspring[i], offspring[j] = method.crossover(p1, p2, rng=rng)
+        else
+            offspring[i], offspring[j] = p1, p2
+        end
+    end
+
+end
+
+function mutate!(population, method, constraints;
+                 rng::AbstractRNG=Random.GLOBAL_RNG)
+    n = length(population)
+    for i in 1:n
+        if rand(rng) < method.mutationRate
+            method.mutation(population[i], rng=rng)
+        end        
+        apply!(constraints, population[i])
+    end
+end
+
+function evaluate!(objfun, fitness, population, constraints)
+    # calculate fitness of the population
+    value!(objfun, fitness, population)
+    # apply penalty to fitness
+    penalty!(fitness, constraints, population)
+end
+

src/mutations.jl

@@ -4,6 +4,13 @@
 # Evolutionary mutations
 # ======================
 
+"""
+    nop(s::AbstractStrategy)
+
+This is a dummy mutation operator that does not change recombinant.
+"""
+nop(recombinant::AbstractVector, s::AbstractStrategy; kwargs...) = recombinant
+
 """
     gaussian(x, s::IsotropicStrategy)
 
@@ -57,6 +64,14 @@ end
 # Strategy mutation operators
 # ===========================
 
+"""
+    nop(s::AbstractStrategy)
+
+This is a dummy operator that does not change strategy.
+"""
+nop(s::AbstractStrategy; kwargs...) = s
+
+
 """
     gaussian(s::IsotropicStrategy)
 
@@ -411,7 +426,7 @@ end
 # ======================
 
 """
-    subtree(method::TreeGP; growth::Real = 0.1)
+    subtree(method::TreeGP; growth::Real = 0.0)
 
 Returns an in-place expression mutation function that performs mutation of an arbitrary expression subtree with a randomly generated one [^5].
 
@@ -523,6 +538,7 @@ function swap!(v::T, from::Int, to::Int) where {T <: AbstractVector}
 end
 
 function mutationwrapper(gamutation::Function)
-    wrapper(recombinant::T, s::S) where {T <: AbstractVector, S <: AbstractStrategy} =  gamutation(recombinant)
+    wrapper(recombinant::T, s::S; kwargs...) where {T <: AbstractVector, S <: AbstractStrategy} =  gamutation(recombinant; kwargs...)
     return wrapper
 end
+

src/nsga2.jl

@@ -67,36 +67,22 @@ function initial_state(method::NSGA2, options, objfun, parents)
 end
 
 function update_state!(objfun, constraints, state, parents::AbstractVector{IT}, method::NSGA2, options, itr) where {IT}
-    @unpack populationSize,crossoverRate,mutationRate,selection,crossover,mutation = method
+    populationSize = method.populationSize
+    rng = options.rng
 
-    # Select offspring
+    # select offspring
     specFit = StackView(state.ranks, state.crowding, dims=1)
-    selected = selection(view(specFit, :, 1:populationSize), populationSize)
-
-    # Perform mating
-    offidx = randperm(populationSize)
-    for i in 1:2:populationSize
-        j = (i == populationSize) ? i-1 : i+1
-        if rand() < crossoverRate
-            state.offspring[i], state.offspring[j] = crossover(parents[selected[offidx[i]]], parents[selected[offidx[j]]])
-        else
-            state.offspring[i], state.offspring[j] = parents[selected[i]], parents[selected[j]]
-        end
-    end
+    selected = method.selection(view(specFit,:,1:populationSize), populationSize; rng=rng)
 
-    # Perform mutation
-    for i in 1:populationSize
-        if rand() < mutationRate
-            mutation(state.offspring[i])
-        end
-        apply!(constraints, state.offspring[i])
-    end
+    # perform mating
+    recombine!(state.offspring, parents, selected, method)
+
+    # perform mutation
+    mutate!(state.offspring, method, constraints, rng=rng)
 
     # calculate fitness of the offspring
     offfit = @view state.fitpop[:, populationSize+1:end]
-    value!(objfun, offfit, state.offspring)
-    # apply penalty to fitness
-    penalty!(offfit, constraints, state.offspring)
+    evaluate!(objfun, offfit, state.offspring, constraints)
 
     # calculate ranks & crowding for population
     F = nondominatedsort!(state.ranks, state.fitpop)
@@ -117,6 +103,7 @@ function update_state!(objfun, constraints, state, parents::AbstractVector{IT},
     state.fittest = state.population[F[1]]
     # and keep their fitness
     state.fitness = state.fitpop[:,F[1]]
+
     # construct new parent population
     parents .= state.population[fitidx]
 

test/knapsack.jl

@@ -1,26 +1,27 @@
 @testset "Knapsack" begin
-    Random.seed!(42)
+    rng = StableRNG(42)
 
     mass    = [1, 5, 3, 7, 2, 10, 5]
     utility = [1, 3, 5, 2, 5,  8, 3]
 
     fitnessFun = n -> (sum(mass .* n) <= 20) ? sum(utility .* n) : 0
 
-    initpop = collect(rand(Bool,length(mass)))
-
+    Random.seed!(rng, 42)
+    initpop = ()->rand(rng, Bool, length(mass))
     result = Evolutionary.optimize(
         x -> -fitnessFun(x),
         initpop,
         GA(
-            selection = tournament(3),
-            mutation = inversion,
+            selection = roulette,
+            mutation = swap2,
             crossover = SPX,
-            mutationRate = 0.9,
-            crossoverRate = 0.2,
-            ɛ = 0.1,                                # Elitism
+            mutationRate = 0.05,
+            crossoverRate = 0.85,
+            ɛ = 0.05,                                # Elitism
             populationSize = 100,
-        ));
-    println("GA:RLT:INV:SP (-objfun) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
+           ), Evolutionary.Options(show_trace=false, rng=rng)
+       );
+    println("GA:RLT:SWP:SPX (-objfun) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
     @test abs(Evolutionary.minimum(result)) == 21.
     @test sum(mass .* Evolutionary.minimizer(result)) <= 20
 
@@ -33,21 +34,22 @@
     uc   = [20]   # upper bound for constraint function
     con = WorstFitnessConstraints(Int[], Int[], lc, uc, cf)
 
-    initpop = BitVector(rand(Bool,length(mass)))
-
+    Random.seed!(rng, 42)
+    initpop = BitVector(rand(rng, Bool, length(mass)))
     result = Evolutionary.optimize(
         x -> -fitnessFun(x), con,
         initpop,
         GA(
-            selection = roulette,
-            mutation = flip,
+            selection = tournament(3),
+            mutation = inversion,
             crossover = SPX,
-            mutationRate = 0.9,
-            crossoverRate = 0.1,
+            mutationRate = 0.05,
+            crossoverRate = 0.85,
             ɛ = 0.1,                                # Elitism
             populationSize = 50,
-        ));
-    println("GA:RLT:FL:SP (-objfun) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
+           ), Evolutionary.Options(show_trace=false, rng=rng, successive_f_tol=10));
+    println("GA:TRN3:INV:SPX (-objfun) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
     @test abs(Evolutionary.minimum(result)) == 21.
     @test sum(mass .* Evolutionary.minimizer(result)) <= 20
+
 end