1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

#light
#r "FSharp.PowerPack.dll"
open System
open Microsoft.FSharp.Math

// interface for test-functions
type TestFunction = 
    abstract Dimension: int 
    abstract Range: (float * float)
    abstract Evaluate: Vector<float> -> float
    
// implementation of rosenbrocks function
type Rosenbrock(n) = 
    interface TestFunction with
        member x.Dimension = n
        member x.Range = (-2.048, 2.048)
        member x.Evaluate xs = 
            xs
            |> Seq.pairwise 
            |> Seq.map 
                (fun (x0,x1) -> 
                    let s = x1 - x0*x0
                    let t = 1.0-x0
                    100.0 * s * s + t*t)
            |> Seq.sum


let rosenbrock100 = new Rosenbrock(100) :> TestFunction


// infinite stream of uniform random numbers 0..1
let uniformRandomStream = 
    let random = new Random()
    Seq.init_infinite (fun i -> random.NextDouble())


// Marsaglia polar method to generate infinite stream of N(0,1) random numbers
let normalRandomStream = 
    let scaledUniformStream = Seq.map (fun x -> 2.0 * x - 1.0 ) uniformRandomStream
    Seq.zip scaledUniformStream scaledUniformStream
    |> Seq.map (fun (x,y) -> (x*x + y*y, x, y))
    |> Seq.filter (fun (s,_,_) -> s < 1.0)
    |> Seq.map_concat (fun (s,x,y) -> 
        let f = sqrt (-2.0 * log s / s)
        [x * f; y * f])


let createRandomSolution (problem:TestFunction) = 
    let min = problem.Range |> (fun (min, _) -> min)
    let max = problem.Range |> (fun (_, max) -> max)
    let range = max - min 
    let solution = 
        uniformRandomStream
        |> Seq.map (fun x -> x*range + min)
        |> Seq.take problem.Dimension
        |> Vector.of_seq
    (solution, 1.0, problem.Evaluate solution)


// adapt sigma and add N(0,s) distributed vector to solution vector
let mutateSolution (problem:TestFunction) tau (xs,s,_) =
    let newS = s * exp (tau * (Seq.hd normalRandomStream))
    let offset = 
        normalRandomStream
        |> Seq.map (fun x -> newS * x)
        |> Seq.take problem.Dimension
        |> Vector.of_seq
    let newXs = xs + offset;
    (newXs, newS, problem.Evaluate newXs)
    

// recombine a seq of solutions by calculating the average vector
let recombineSolutions selected = 
    let n = float (Seq.length selected)
    let (newX,newS,_) = 
        let x0 = Seq.hd selected
        Seq.skip 1 selected
        |> Seq.fold 
            (fun (x:vector,sx,_) (y:vector,sy,_) -> (x+y,sx+sy,Double.NaN)) x0
    (newX * (1.0/n), newS/n, Double.NaN)
   

let evolutionStrategy (problem:TestFunction) iterations mu rho lambda = 
    let tau = 1.0 / sqrt (float problem.Dimension) // learning rate (for sigma adaption)
    let sort pop = 
        Array.sort_by (fun (_,_,quality) -> quality) pop
        pop
    let select lambda rho pop =                    // select solutions from pop randomly
        let max = float mu
        let selectRecombinants rho pop =
            uniformRandomStream 
            |> Seq.map (fun u -> Array.get pop (int (floor (u * max))))
            |> Seq.take rho            
        Seq.init_finite lambda (fun i -> (selectRecombinants rho pop))
    let recombine = Seq.map recombineSolutions
    let mutate = Seq.map (fun x -> (mutateSolution problem tau x))
    let replace pop children =                     // keep mu best from offspring and parents (mu+lambda)
        Seq.append pop children
        |> Array.of_seq
        |> sort
        |> (fun arr -> Array.sub arr 0 mu)
    let printStats pop = 
        pop |> Seq.hd |> (fun (_,_,quality) -> printf "Quality %f\n" quality) 
        pop
    let rec iterate env step n = 
      if n<=0 then env else iterate (step env) step (n - 1)
    iterate
        ([|for i in 1..mu -> createRandomSolution problem|] |> sort) // init random population
        (fun pop ->
            pop |> printStats |> select lambda rho |> recombine |> mutate |> replace pop) // main loop
        iterations


// find solution for 100-dimensional rosenbrock function
// with ES(100/2, 1000) 
// 100 iterations
evolutionStrategy  rosenbrock100 100 100 2 1000 |> Seq.hd