Differential evolution in Go

evo.go

// Keep updated here https://gist.github.com/peheje/ea1d067803c8bb2a9fdbbc17b763174a

package main

import (
	"fmt"
	"math"
	"math/rand"
	"time"

	"golang.org/x/exp/constraints"
)

func f1(xs []float64) float64 {
	s := 0.0
	for i := 0; i < len(xs); i++ {
		s += (xs[i] * xs[i])
	}
	return s
}

func f2(xs []float64) float64 {
	s := 0.0
	p := 1.0
	for i := 0; i < len(xs); i++ {
		s += math.Abs(xs[i])
		p *= math.Abs(xs[i])
	}
	return math.Abs(s) + math.Abs(p)
}

func booth(xs []float64) float64 {
	t1 := xs[0] + 2.0*xs[1] - 7.0
	t2 := 2.0*xs[0] + xs[1] - 5.0
	return t1*t1 + t2*t2
}

func nonlinearSystemOfEquationsWithConstraints(xs []float64) float64 {
	x := xs[0]
	y := xs[1]
	t1 := (x + 1.0) * (10.0 - x) * ((1.0 + y*y) / (1.0 + y*y + y))
	t2 := (y + 2.0) * (20.0 - y) * ((1.0 + x*x) / (1.0 + x*x + x))
	c1 := x*x + y*y

	e := 0.0
	if c1 > 10.0 {
		e = c1
	}

	return math.Abs(t1) + math.Abs(t2) + e
}

func rastrigin(xs []float64) float64 {
	n := float64(len(xs))
	a := 10.0

	s := 0.0
	for _, x := range xs {
		s += x*x - (a * math.Cos(2.0*math.Pi*x))
	}

	return a*n + s
}

func randf(min, max float64) float64 {
	return min + rand.Float64()*(max-min)
}

func randi(max int) int {
	return rand.Intn(max)
}

func clamp(value, min, max float64) float64 {
	if value > max {
		return max
	} else if value < min {
		return min
	} else {
		return value
	}
}

type Number interface {
	constraints.Float | constraints.Integer
}

func average[T Number](xs []T) float64 {
	s := 0.0
	for _, x := range xs {
		s += float64(x)
	}
	return s / float64(len(xs))
}

func minimum[T Number](xs []T) (int, T) {
	min := xs[0]
	mini := 0
	for i := 1; i < len(xs); i++ {
		if xs[i] < min {
			min = xs[i]
			mini = i
		}
	}

	return mini, min
}

func main() {

	start := time.Now()
	print := 1000
	optimizer := rastrigin
	argsize := 200
	boundsMin := -10.0
	boundsMax := 10.0
	generations := 10_000
	popsize := 200
	mutateMin := 0.2
	mutateMax := 0.95
	crossoverMin := 0.1
	crossoverMax := 1.0

	trial := make([]float64, argsize)
	pop := make([][]float64, popsize)

	for i := range pop {
		pop[i] = make([]float64, argsize)
		for j := range pop[i] {
			pop[i][j] = randf(boundsMin, boundsMax)
		}
	}

	scores := make([]float64, popsize)
	for i := range scores {
		scores[i] = f1(pop[i])
	}

	for g := 0; g < generations; g++ {
		crossover := randf(crossoverMin, crossoverMax)
		mutate := randf(mutateMin, mutateMax)

		for i := 0; i < popsize; i++ {
			x0 := pop[randi(popsize-1)]
			x1 := pop[randi(popsize-1)]
			x2 := pop[randi(popsize-1)]
			xt := pop[i]

			for j := 0; j < argsize; j++ {
				if randf(0.0, 1.0) < crossover {
					trial[j] = clamp((x0[j]+(x1[j]-x2[j]))*mutate, boundsMin, boundsMax)
				} else {
					trial[j] = xt[j]
				}
			}

			scoreTrial := optimizer(trial)
			if scoreTrial < scores[i] {
				copy(pop[i], trial)
				scores[i] = scoreTrial
			}
		}

		if g%print == 0 || g == generations-1 {
			_, min := minimum(scores)
			mean := average(scores)

			fmt.Printf("generation %d\n", g)
			fmt.Printf("generation mean %v\n", mean)
			fmt.Printf("generation best %v\n", min)
		}
	}

	minIndex, min := minimum(scores)
	fmt.Printf("best %v\n", pop[minIndex])
	fmt.Printf("generation best %v\n", min)
	end := time.Now()
	fmt.Printf("time taken: %vs\n", (end.Sub(start).Seconds()))
}