allow better control over used RNG (part 1)

- added `rng` parameter to `Options`
- wired calls that use RNG to use the option `rng` parameter in GA & GP

Author: wildart

src/api/types.jl

@@ -36,6 +36,7 @@ There are following options available:
     callback::TCallback = nothing
     time_limit::Float64 = NaN
     parallelization::Symbol = :serial
+    rng::AbstractRNG = Random.GLOBAL_RNG
 end
 function show(io::IO, o::Options)
     for k in fieldnames(typeof(o))

src/ga.jl

@@ -61,19 +61,19 @@ end
 function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::GA, options, itr) where {IT}
     @unpack populationSize,crossoverRate,mutationRate,ɛ,selection,crossover,mutation = method
     evaltype = options.parallelization
-
+    rng = options.rng
     offspring = similar(population)
 
     # Select offspring
-    selected = selection(state.fitpop, populationSize)
+    selected = selection(state.fitpop, populationSize, rng=rng)
 
     # Perform mating
-    offidx = randperm(populationSize)
+    offidx = randperm(rng, populationSize)
     offspringSize = populationSize - state.eliteSize
     for i in 1:2:offspringSize
         j = (i == offspringSize) ? i-1 : i+1
-        if rand() < crossoverRate
-            offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]])
+        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
@@ -88,8 +88,8 @@ function update_state!(objfun, constraints, state, population::AbstractVector{IT
 
     # Perform mutation
     for i in 1:offspringSize
-        if rand() < mutationRate
-            mutation(offspring[i])
+        if rand(rng) < mutationRate
+            mutation(offspring[i], rng=rng)
         end
     end
 

src/gp.jl

@@ -55,7 +55,7 @@ function randterm(rng::AbstractRNG, t::TreeGP)
     if isa(term, Symbol) || isa(term, Real)
         term
     elseif isa(term, Function)
-        term()
+        term(rng) # terminal functions must accept RNG as an argument
     else
         # Code shouldn't reach branch but left as a catchall
         dim = t.terminals[term]
@@ -97,10 +97,11 @@ rand(t::TreeGP, maxdepth::Int=2; kwargs...) =
 
 Initialize a random population of expressions derived from `expr`.
 """
-function initial_population(m::TreeGP, expr::Union{Expr,Nothing}=nothing)
+function initial_population(m::TreeGP, expr::Union{Expr,Nothing}=nothing;
+                            rng::AbstractRNG=Random.GLOBAL_RNG)
     n = population_size(m)
     return [
-        expr === nothing ? rand(m, m.maxdepth, mindepth=m.mindepth) : deepcopy(expr)
+        expr === nothing ? rand(rng, m, m.maxdepth, mindepth=m.mindepth) : deepcopy(expr)
         for i in 1:n
     ]
 end
@@ -129,8 +130,13 @@ function update_state!(objfun, constraints, state::GPState, population::Abstract
 end
 
 # Custom optimization call
-function optimize(f, method::TreeGP, options::Options = Options(;default_options(method)...))
-    method.optimizer.mutation =  subtree(method)
-    method.optimizer.crossover = crosstree
-    optimize(f, NoConstraints(), nothing, method, options)
+function optimize(f, mthd::TreeGP, options::Options = Options(;default_options(mthd)...))
+    optimize(f, mthd, initial_population(mthd, rng=options.rng), options)
 end
+
+function optimize(f, mthd::TreeGP, population, options::Options = Options(;default_options(mthd)...))
+    mthd.optimizer.mutation =  subtree(mthd)
+    mthd.optimizer.crossover = crosstree
+    optimize(f, NoConstraints(), mthd, population, options)
+end
+

test/gp.jl

@@ -1,5 +1,6 @@
 @testset "Genetic Programming" begin
-    Random.seed!(9874984737486);
+    rng = StableRNG(42)
+
     pop = 10
     terms = Terminal[:x, :y, rand]
     funcs = Function[+,-,*,/]
@@ -16,18 +17,18 @@
     @test_skip summary(t) == "TreeGP[P=10,Parameter[x,y],Function[*, +, /, -]]"
 
     # population initialization
-    popexp = Evolutionary.initial_population(t);
+    popexp = Evolutionary.initial_population(t, rng=rng);
     @test length(popexp) == pop
     popexp = Evolutionary.initial_population(t, :(x + 1));
     @test popexp[1] == :(x + 1)
 
     # recursive helper functions
-    Random.seed!(9874984737482)
-    gt = rand(TreeGP(pop, terms, funcs, maxdepth=2, initialization=:grow), 3)
+    Random.seed!(rng, 1)
+    gt = rand(rng, TreeGP(pop, terms, funcs, maxdepth=2, initialization=:grow), 3)
     @test Evolutionary.nodes(gt) < 15
     @test Evolutionary.height(gt) <= 3
     @test length(gt) < 15
-    ft = rand(TreeGP(pop, terms, funcs, maxdepth=2, initialization=:full), 3)
+    ft = rand(rng, TreeGP(pop, terms, funcs, maxdepth=2, initialization=:full), 3)
     @test Evolutionary.nodes(ft) == 15
     @test Evolutionary.height(ft) == 3
     @test length(ft) == 15
@@ -57,7 +58,7 @@
 
     # evaluation
     ex = Expr(:call, +, 1, :x) |> Evolutionary.Expression
-    xs = rand(10)
+    xs = rand(rng, 10)
     @test ex(xs[1]) == xs[1]+1
     @test ex.(xs) == xs.+1
     io = IOBuffer()
@@ -67,24 +68,23 @@
     depth = 5
     fitfun(x) = x*x + x + 1.0
     function fitobj(expr)
-        rg = -5.0:0.5:5.0
+        rg = -5.0:0.1:5.0
         ex = Evolutionary.Expression(expr)
-        sum(v->isnan(v) ? 1.0 : v, abs2.(fitfun.(rg) - ex.(rg)) )/length(rg) |> sqrt
+        sum(v->isnan(v) ? 1.0 : v, abs2.(fitfun.(rg) - ex.(rg)) )/length(rg)
     end
 
-    Random.seed!(9874984737426);
     res = Evolutionary.optimize(fitobj,
         TreeGP(50, Terminal[:x, randn], Function[+,-,*,Evolutionary.pdiv],
             mindepth=1,
             maxdepth=depth,
             optimizer = GA(
-                selection = tournament(5),
+                selection = tournament(4),
                 ɛ = 0.1,
-                mutationRate = 0.85,
+                mutationRate = 0.8,
                 crossoverRate = 0.2,
             ),
-        )
+        ), Evolutionary.Options(rng=StableRNG(1))
     )
-    @test minimum(res) < 1.1
+    @test minimum(res) < 1
 
 end

test/rosenbrock.jl

@@ -52,28 +52,28 @@
     c = PenaltyConstraints(100.0, fill(0.0, 5N), Float64[], [1.0], [1.0], con_c!)
     result = Evolutionary.optimize(rosenbrock, c, (() -> rand(5N)), CMAES(mu = 40, lambda = 100))
     println("(5/5,100)-CMA-ES [penalty] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
-    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=1e-1
+    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=0.1
 
     c = PenaltyConstraints(100.0, fill(0.0, 2N), fill(0.5, 2N), [1.0], [1.0], con_c!)
     result = Evolutionary.optimize(rosenbrock, c, (() -> rand(2N)), CMAES(mu = 8, lambda = 100))
     println("(5/5,100)-CMA-ES [penalty] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
-    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=1e-1
+    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=0.1
     @test all(0.0 <= x+0.01 && x-0.01 <= 0.5 for x in abs.(Evolutionary.minimizer(result)))
 
     # Testing: GA
     m = GA(
         populationSize = 100,
         ɛ = 0.1,
         selection = rouletteinv,
-        crossover = intermediate(0.25),
+        crossover = IC(0.2),
         mutation = BGA(fill(0.5,N))
     )
     result = Evolutionary.optimize(rosenbrock, (() -> rand(N)), m)
     println("GA(p=100,x=0.8,μ=0.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
-    test_result(result, N, 1e-1)
+    test_result(result, N, 0.1)
     result = Evolutionary.optimize(rosenbrock, BoxConstraints(0.0, 0.5, N), (() -> rand(N)), m)
     println("GA(p=100,x=0.8,μ=0.1,ɛ=0.1)[box] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
-    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=1e-1
+    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=0.1
 
     # Testing: DE
     result = Evolutionary.optimize(rosenbrock, (() -> rand(N)), DE(populationSize = 100))