added new combinatorial stuff:

- crossover: subset selection (SSX)
- mutation: `replace`
- problem: 120 bits subset selection

Author: wildart

docs/src/crossover.md

@@ -10,6 +10,8 @@ All recombination operations have following call interface: `recombination(i1, i
 
 ## Operations
 
+### ES Crossovers
+
 List of the ES strategy recombination operations:
 
 ```@docs
@@ -23,7 +25,7 @@ average(population::Vector{T}) where {T <: AbstractVector}
 marriage
 ```
 
-Binary crossovers:
+### Binary Crossovers
 
 ```@docs
 SPX
@@ -35,7 +37,7 @@ EXPX
 BSX
 ```
 
-Real-valued crossovers:
+### Real-valued Crossovers
 
 ```@docs
 identity
@@ -50,17 +52,18 @@ MILX
 SBX
 ```
 
-Combinatorial crossovers:
+### Combinatorial Crossovers
 
 ```@docs
 PMX
 OX1
 CX
 OX2
 POS
+SSX
 ```
 
-Tree (expression) crossovers:
+### Tree (expression) Crossovers
 
 ```@docs
 Evolutionary.crosstree
@@ -80,4 +83,4 @@ Evolutionary.crosstree
 
 [^6]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
 
-[^7]: M. A. Wolters, “A Genetic Algorithm for Selection of Fixed-Size Subsets with Application to Design Problems”, J. Stat. Soft., vol. 68, no. 1, pp. 1–18, Nov. 2015.
+[^7]: M. A. Wolters, “A Genetic Algorithm for Selection of Fixed-Size Subsets with Application to Design Problems”, J. Stat. Soft., vol. 68, no. 1, pp. 1–18, Nov. 2015.

docs/src/mutation.md

@@ -57,6 +57,7 @@ insertion
 swap2
 scramble
 shifting
+Evolutionary.replace
 ```
 
 ## Differential Evolution
@@ -92,4 +93,4 @@ shrink
 
 [^8]: P. J. Angeline, "An investigation into the sensitivity of genetic programming to the frequency of leaf selection during subtree crossover", Genetic Programming 1996: Proceedings of the First Annual Conference, 21–29, 1996.
 
-[^9]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
+[^9]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995

examples/120bits.jl

@@ -0,0 +1,76 @@
+"""
+Problem: 120-bit sparse subset
+
+Let a 120-bit string be b = [b_1|b_2|···|b_20], where b_i, i = 1, 2, ..., 20 are
+six-bit substrings referred to as blocks. Define the sequence `000001` to be
+the target configuration for a block, and denote this configuration by b_t.
+The objective function is defined such that the optimal solution is
+b* = [b_t|b_t|···|b_t], that is, the string with all 20 blocks having the target
+configuration.
+
+This solution can be expressed as the index vector v* = [6, 12, 18,...,120].
+"""
+
+using Evolutionary, Test, Random
+
+
+function fitness(oset, bsize = 6)
+    scores = Dict(i => (i != bsize ? 3i : 20) for i in 0:bsize)
+    # count 1s in blocks
+    ht = Dict{Int,Int}()
+    sl = sr = 0
+    for i in oset
+        if i > 60
+            sr += 1
+        else
+            sl += 1
+        end
+        blk = div(i, bsize)
+        blk -= i % bsize == 0 ? 1 : 0
+        if haskey(ht, blk)
+            ht[blk] += 1
+        else
+            ht[blk] = 1
+        end
+    end
+    # calculate total score
+    ss = 0
+    for blk in sort(collect(keys(ht)))
+        cnt = ht[blk]
+        if cnt == 7
+            println(oset)
+            println(ht)
+            error("Problem")
+        end
+        ss += scores[cnt]
+        if (cnt == 1) && ((blk+1)*bsize in oset)
+            ss += 9
+        end
+        #println("$blk, $cnt, $(scores[cnt]), $ss")
+    end
+    ss -= 2(sl > 13 ? sl : 0)
+    ss -= 2(sr > 13 ? sr : 0)
+    ss
+end
+
+# Bit string examples 
+bstr = "000100 001010 000010 111111 000100 101000 000010 000000 000001 000000 000010 000001 001000 000000 000000 010000 000000 000010 000000 000000"
+idxs = [i for (i,j) in enumerate(filter(c->c != ' ', bstr)) if j == '1']
+@test fitness(idxs) == 50
+
+# optimal solution
+idxs_opt = [6*i for i in 1:20]
+@test fitness(idxs_opt) == 240
+
+# GA solution
+mthd = GA(populationSize = 100, selection = tournament(5),
+          crossover = SSX, crossoverRate = 0.99,
+          mutation  = replace(collect(1:120)), mutationRate = 0.1)
+
+Random.seed!(42)
+opts = Evolutionary.Options(show_trace=false, successive_f_tol=30)
+init_pop = ()->randperm(120)[1:20]
+res = Evolutionary.optimize(idxs->-fitness(idxs), init_pop, mthd, opts)
+println(res)
+Evolutionary.minimizer(res) |> sort
+

src/Evolutionary.jl

@@ -6,7 +6,7 @@ module Evolutionary
     using NLSolversBase: NLSolversBase, AbstractObjective, ConstraintBounds,
                          AbstractConstraints, nconstraints_x, nconstraints
     import NLSolversBase: f_calls, value, value!, value!!
-    import Base: show, copy, minimum, summary, identity, getproperty, rand, getindex, length, copyto!, setindex!
+    import Base: show, copy, minimum, summary, identity, getproperty, rand, getindex, length, copyto!, setindex!, replace
 
     export AbstractStrategy, strategy, mutationwrapper,
            IsotropicStrategy, AnisotropicStrategy, NoStrategy,
@@ -17,14 +17,14 @@ module Evolutionary
            gaussian, cauchy,
            # GA mutations
            flip, bitinversion, uniform, BGA, inversion, insertion, swap2, scramble,
-           shifting, PM, MIPM, PLM,
+           shifting, PM, MIPM, PLM, replace,
            # ES recombinations
            average, marriage,
            # GA recombinations
            SPX, TPX, UX, SHFX,
            DC, AX, WAX, IC, LC, HX, LX, MILX, SBX,
            PMX, OX1, CX, OX2, POS,
-           BSX,
+           BSX, SSX,
            # GA selections
            ranklinear, uniformranking, roulette, rouletteinv, sus, susinv,
            tournament, truncation,

src/mutations.jl

@@ -356,6 +356,36 @@ function shifting(recombinant::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <
     return recombinant
 end
 
+"""
+    replace(pool,[minchange=1])(recombinant)
+
+Replacement mutation operator changes an arbitrary number, no smaller then `minchange`,
+of elements in the individual by replacing them with elements from the predefined `pool` that are not in the individual.
+"""
+function replace(pool::Vector{P}; minchange=1) where {P}
+    function rplc(recombinant::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <: AbstractVector}
+        l = length(recombinant)
+        p = length(pool)
+        # how many values to change
+        nchg = max(minchange, rand(rng, 1:min(l, p-l)))
+        # select new elements
+        idxs = randperm(rng, p)
+        new_vals = P[]
+        for i in idxs
+            if pool[i] ∉ recombinant
+                push!(new_vals, pool[i])
+            end
+            length(new_vals) == nchg && break
+        end
+        # update arbitrary positions with new values
+        new_idxs = randperm(rng, l)[1:nchg]
+        recombinant[new_idxs] .= new_vals
+        return recombinant
+    end
+    return rplc
+end
+
+
 # Differential Evolution
 # ======================
 
@@ -376,6 +406,7 @@ function differentiation(recombinant::T, mutators::AbstractVector{T}; F::Real =
     return recombinant
 end
 
+
 # Genetic Programming
 # ======================
 

src/recombinations.jl

@@ -95,7 +95,7 @@ end
 """
     SHFX(v1, v2)
 
-Shuffle crossover between the parents `v1` and `v2` that performs recombination similar to (SPX)[@ref] preliminary shuffling these parents.
+Shuffle crossover between the parents `v1` and `v2` that performs recombination similar to [`SPX`](@ref) preliminary shuffling these parents.
 """
 function SHFX(v1::T, v2::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <: AbstractVector}
     l = length(v1)
@@ -132,7 +132,7 @@ end
 """
     BINX(Cr::Real=0.5)
 
-Returns a uniform (binomial) crossover function, see [Recombination Interface](@ref), function with the propbabilty `Cr` [^2].
+Returns a uniform (binomial) crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].
 
 The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
 """
@@ -155,7 +155,7 @@ end
 """
     EXPX(Cr::Real=0.5)
 
-Returns an exponential crossover function, see [Recombination Interface](@ref), function with the propbabilty `Cr` [^2].
+Returns an exponential crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].
 
 The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
 """
@@ -183,7 +183,7 @@ end
 """
     BSX(k::Int)
 
-Binary Subset Crossover[^7]. Produces two offsprings by first pooling the unique
+Binary Subset Crossover[^7]. Produces an offspring by first pooling the unique
 items of the two parents, and then creating each offspring by sampling without 
 replacement at most `k` elements from the pool of items.
 """
@@ -201,6 +201,7 @@ function BSX(k::Int)
     return BSX
 end
 
+
 # Real valued crossovers
 # ----------------------
 
@@ -533,7 +534,7 @@ end
 """
     POS(v1, v2)
 
-Position-based crossover is a modification of the [OX1](@ref) operator.
+Position-based crossover is a modification of the [`OX1`](@ref) operator.
 It selects a random set of positions in the parents `v1` and `v2`, then imposes
 the position of the selected elements of one parent on the corresponding elements
 of the other parent.
@@ -549,7 +550,7 @@ function POS(v1::T, v2::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <: Abstr
             c1[i] = v2[i]
             c2[i] = v1[i]
         end
-    end
+end
 
     for i in 1:s
         if !in(v1[i],c1)
@@ -564,6 +565,24 @@ function POS(v1::T, v2::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <: Abstr
     return c1,c2
 end
 
+"""
+    SSX(v1, v2)
+
+Subset crossover operator creates new offspring by pooling the unique indices of
+the two parent vectors `v1` and `v2`, and then sampling a set of unique indices
+from this pool, uniformly at random.
+"""
+function SSX(v1::T, v2::T; rng::AbstractRNG=Random.GLOBAL_RNG) where {T <: AbstractVector}
+    l = length(v1)
+    pool = unique([v1;v2])
+    p = length(pool)
+    c1, c2 = if l <= p
+        pool[randperm(p)[1:l]], pool[randperm(p)[1:l]] 
+    else
+        rand(pool, l), rand(pool, l)
+    end 
+    return c1, c2
+end
 
 
 # ===================

test/gp.jl

@@ -80,7 +80,7 @@
             optimizer = GA(
                 selection = tournament(5),
                 ɛ = 0.1,
-                mutationRate = 0.8,
+                mutationRate = 0.85,
                 crossoverRate = 0.2,
             ),
         )

test/mutations.jl

@@ -62,6 +62,8 @@
         @test scramble([0:5;], rng=rng) == [0,4,1,5,3,2]
         Random.seed!(rng, 2)
         @test shifting([0:5;], rng=rng) == [2,3,4,5,0,1]
+        Random.seed!(rng, 2)
+        @test replace([0:9;])([0:5;], rng=rng) == [7,1,2,3,4,9]
 
     end
 

test/recombinations.jl

@@ -59,6 +59,12 @@
             @test off1[end] == 0 && off2[end] == 0
             @test sum(off1) == 3 && sum(off2) == 3
         end
+
+        Random.seed!(rng, 1)
+        off1, off2 = SSX([0:4;], [5:9;], rng=rng)
+        @test all(i -> i in 0:9, off1)
+        @test all(i -> i in 0:9, off2)
+
     end
 
     @testset "DE" begin