ga: implement SOMA T3A algorithm

algo/ga/error.go

@@ -0,0 +1,28 @@
+// Copyright 2023 wanderer <a_mirre at utb dot cz>
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+package ga
+
+type gaError struct {
+	s string
+}
+
+var (
+	errGenerationsIsZero    = newGAError("Generations cannot be 0")
+	errBenchMinItersIsLTOne = newGAError("Minimum bench iterations cannot be less than 1")
+	errNPIsZero             = newGAError("NP cannot be 0")
+	errKIsZero              = newGAError("K cannot be 0")
+	errMIsZero              = newGAError("M cannot be 0")
+	errNIsZero              = newGAError("N cannot be 0")
+	errNjumpsIsZero         = newGAError("Njumps cannot be 0")
+	errDimensionsUnset      = newGAError("Dimensions cannot be unset")
+	errBenchUnset           = newGAError("Bench cannot be unset")
+)
+
+func (e *gaError) Error() string {
+	return e.s
+}
+
+func newGAError(text string) error {
+	return &gaError{s: text}
+}

algo/ga/somat3a.go

@@ -4,7 +4,12 @@
 package ga
 
 import (
+	"sort"
+	"time"
+
+	"git.dotya.ml/wanderer/math-optim/bench/cec2020"
 	"git.dotya.ml/wanderer/math-optim/stats"
+	"golang.org/x/exp/rand"
 	"gonum.org/v1/gonum/stat/distuv"
 )
 
@@ -21,6 +26,17 @@ type SOMAT3A struct {
 	Dimensions []int
 	// NP is the initial population size.
 	NP int
+	// K is the number of individuals to choose the leader from.
+	K int
+	// M denotes how many individuals are picked from the population to form a
+	// `team` during the organisation phase, can be thought of as "team size".
+	M int
+	// N is the number of best individuals in each team selected for actual
+	// migration.
+	N int
+	// Njumps is the fixed number of jumps that each chosen migrating
+	// individual performs on their way to the leader.
+	Njumps int
 	// BenchName is the human-friendly name of the benchmarking function.
 	BenchName string
 	// ch is a channel for writing back computed results.
@@ -34,3 +50,230 @@ type SOMAT3A struct {
 	// initialised denotes the initialisation state of the struct.
 	initialised bool
 }
+
+var somat3aLogger = newLogger(" ***  SOMA T3A")
+
+// Init performs checks on the input params and initialises an instance of
+// SOMAT3A.
+func (s *SOMAT3A) Init(
+	generations, benchMinIters, np, k, m, n, njumps int,
+	dimensions []int,
+	bench string,
+	ch chan []stats.Stats,
+	chAlgoMeans chan *stats.AlgoBenchMean,
+) error {
+	switch {
+	case generations == 0:
+		return errGenerationsIsZero
+
+	case generations == -1:
+		somat3aLogger.Println("Generations is '-1', disabling generation limits")
+
+	case benchMinIters < 1:
+		return errBenchMinItersIsLTOne
+
+	case np == 0:
+		return errNPIsZero
+
+	case k == 0:
+		return errKIsZero
+
+	case m == 0:
+		return errMIsZero
+
+	case n == 0:
+		return errNIsZero
+
+	case njumps == 0:
+		return errNjumpsIsZero
+
+	case len(dimensions) == 0:
+		return errDimensionsUnset
+
+	case bench == "":
+		return errBenchUnset
+	}
+
+	s.Generations = generations
+	s.BenchMinIters = benchMinIters
+	s.NP = np
+	s.K = k
+	s.M = m
+	s.N = n
+	s.Njumps = njumps
+
+	s.Dimensions = dimensions
+	s.BenchName = bench
+	s.ch = ch
+	s.chAlgoMeans = chAlgoMeans
+
+	s.initialised = true
+
+	gaLogger.Printf("%s init done for bench %s", s.getAlgoName(), s.BenchName)
+
+	return nil
+}
+
+// Run runs the SOMA T3A algo for each of `s.Dimensions` with benchmarking
+// function `s.BenchName`, repeats the run `s.BenchMinIters` times and returns
+// the stats via chans.
+func (s *SOMAT3A) Run() {
+	if !s.initialised {
+		somat3aLogger.Fatalln(s.getAlgoName(), "needs to be initialised before calling Run(), exiting..")
+	}
+
+	var somat3aStats []stats.Stats
+	// testing purposes.
+	// s.BenchMinIters = 5
+
+	somat3aMeans := &stats.AlgoBenchMean{
+		Algo:       s.getAlgoName(),
+		BenchMeans: make([]stats.BenchMean, 0, len(s.Dimensions)),
+	}
+
+	for _, dim := range s.Dimensions {
+		maxFES := cec2020.GetMaxFES(dim)
+		// TODO(me): check the slides for the formula to calculate this.
+		// 4 fes per member per njump.
+		// fesPerIter := int(float64(maxFES / (4 * s.NP * s.Njumps)))
+		fesPerIter := int(float64(maxFES / (4 * s.NP)))
+		somat3aStatDimX := &stats.Stats{
+			Algo:        s.getAlgoName(),
+			Dimens:      dim,
+			Iterations:  s.BenchMinIters,
+			Generations: maxFES,
+			NP:          s.NP,
+		}
+		funcStats := &stats.FuncStats{BenchName: s.BenchName}
+		dimXMean := &stats.BenchMean{
+			Bench:       s.BenchName,
+			Dimens:      dim,
+			Iterations:  s.BenchMinIters,
+			Generations: maxFES,
+			Neighbours:  -1,
+		}
+		uniDist := distuv.Uniform{
+			Min: cec2020.SearchRange.Min(),
+			Max: cec2020.SearchRange.Max(),
+		}
+
+		if maxFES == -1 {
+			somat3aLogger.Fatalf("could not get maxFES for current dim (%d), bailing", dim)
+		}
+
+		somat3aLogger.Printf("running bench \"%s\" for %dD, maxFES: %d\n",
+			s.BenchName, dim, maxFES,
+		)
+
+		funcStats.BenchResults = make([]stats.BenchRound, s.BenchMinIters)
+
+		// rand.Seed(uint64(time.Now().UnixNano()))
+		// create and seed a source of pseudo-randomness
+		src := rand.NewSource(uint64(rand.Int63()))
+
+		// track execution duration.
+		start := time.Now()
+
+		for iter := 0; iter < s.BenchMinIters; iter++ {
+			somat3aLogger.Printf("run: %d, bench: %s, %dD, started at %s",
+				iter, s.BenchName, dim, start,
+			)
+
+			funcStats.BenchResults[iter].Iteration = iter
+
+			uniDist.Src = src
+			// uniDist.Src = rand.NewSource(uint64(rand.Int63()))
+
+			// create a population with known params.
+			pop := newSOMAT3APopulation(s.BenchName, s.NP, s.K, s.M, s.N, s.Njumps, dim)
+
+			// set population seed.
+			pop.Seed = uint64(time.Now().UnixNano())
+			// initialise the population.
+			pop.Init()
+
+			bestResult := s.evaluate(pop)
+
+			// the core.
+			for i := 0; i < fesPerIter; i++ {
+				funcStats.BenchResults[iter].Results = append(
+					funcStats.BenchResults[iter].Results,
+					bestResult,
+				)
+
+				// call evolve where the inner workngs of SOMA T3A run.
+				s.evolve(pop)
+
+				// take measurements of current population fitness.
+				r := s.evaluate(pop)
+
+				// save if better.
+				if r < bestResult {
+					bestResult = r
+				}
+
+				// this block makes sure we properly count func evaluations for
+				// the purpose of correctly comparable plot comparison. i.e.
+				// append the winning (current best) value 4*NP-1 (the first
+				// best is already saved at this point) times to represent the
+				// fact that while evaluating (and comparing) other population
+				// individuals to the current best value is taking place in the
+				// background, the current best value itself is kept around and
+				// symbolically saved as the best of the Generation.
+				for x := 0; x < 4*s.NP-1; x++ {
+					funcStats.BenchResults[iter].Results = append(
+						funcStats.BenchResults[iter].Results,
+						bestResult,
+					)
+				}
+			}
+		}
+
+		elapsed := time.Since(start)
+
+		somat3aLogger.Printf("completed: bench: %s, %dD, computing took %s\n",
+			s.BenchName, dim, elapsed,
+		)
+
+		// get mean vals.
+		dimXMean.MeanVals = stats.GetMeanVals(funcStats.BenchResults, maxFES)
+		funcStats.MeanVals = dimXMean.MeanVals
+
+		somat3aStatDimX.BenchFuncStats = append(
+			somat3aStatDimX.BenchFuncStats, *funcStats,
+		)
+
+		somat3aStats = append(somat3aStats, *somat3aStatDimX)
+
+		somat3aMeans.BenchMeans = append(somat3aMeans.BenchMeans, *dimXMean)
+	}
+
+	sort.Sort(somat3aMeans)
+
+	s.chAlgoMeans <- somat3aMeans
+	s.ch <- somat3aStats
+}
+
+func (s *SOMAT3A) evaluate(pop *SOMAT3APopulation) float64 {
+	bestIndividual := pop.Population[pop.GetBestIdx()]
+
+	return pop.ProblemFunc(bestIndividual.CurX)
+}
+
+func (s *SOMAT3A) evolve(pop *SOMAT3APopulation) {
+	// organise.
+	leader, migrants := pop.organise()
+	// migrate.
+	pop.migrate(leader, migrants)
+	// actualise (position, best values..).
+	pop.actualise(migrants)
+}
+
+func (s *SOMAT3A) getAlgoName() string {
+	return "SOMA T3A"
+}
+
+// NewSOMAT3A returns a pointer to a new, uninitialised SOMAT3A instance.
+func NewSOMAT3A() *SOMAT3A {
+	return &SOMAT3A{}
+}

algo/ga/somat3aActualisation.go

@@ -0,0 +1,31 @@
+// Copyright 2023 wanderer <a_mirre at utb dot cz>
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+package ga
+
+func (p *SOMAT3APopulation) actualise(migrants []int) {
+	f := p.ProblemFunc
+
+	for _, i := range migrants {
+		bestJump := p.Population[i].Jumps[0]
+
+		for j := 1; j < p.Njumps; j++ {
+			jump := p.Population[i].Jumps[j]
+
+			// incrementing fes for the following is done in
+			// somat3aMigration.go.
+			if f(jump) < f(bestJump) {
+				bestJump = jump
+			}
+		}
+
+		if f(bestJump) < f(p.Population[i].CurX) {
+			p.Population[i].CurX = bestJump
+		}
+
+		p.fes += 2
+
+		// clean jump history
+		p.Population[i].Jumps = nil
+	}
+}

algo/ga/somat3aMigration.go

@@ -0,0 +1,80 @@
+// Copyright 2023 wanderer <a_mirre at utb dot cz>
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+package ga
+
+import (
+	"git.dotya.ml/wanderer/math-optim/bench/cec2020"
+)
+
+// migration is used in SOMA-type algos.
+
+func (p *SOMAT3APopulation) migrate(leader int, migrants []int) {
+	for _, i := range migrants {
+		jumps := make([][]float64, p.Njumps)
+
+		// jump zero.
+		// jumps[0] = p.jump(i, leader, getStep(0, maxFES))
+
+		// bestJump of this migrant is now the jump zero.
+		// bestJump := jumps[0]
+
+		for j := 0; j < p.Njumps; j++ {
+			step := p.getStep()
+			jump := p.jump(i, leader, step)
+
+			p.fes += 2
+
+			jumps[j] = jump
+		}
+
+		p.Population[i].Jumps = jumps
+
+		p.Population[i].PRTVector = p.getPrtVector()
+	}
+
+	p.PRT = p.getPrt()
+}
+
+func (p *SOMAT3APopulation) jump(migrant, leader int, step float64) []float64 {
+	m := p.Population[migrant]
+	l := p.Population[leader]
+	prtv := m.PRTVector
+	t := step
+	jumpPosition := make([]float64, len(m.CurX))
+
+	for i := range m.CurX {
+		jumpPosition[i] = m.CurX[i] + ((l.CurX[i] - m.CurX[i]) * t * prtv[i])
+	}
+
+	return jumpPosition
+}
+
+// getPrtVector creates and returns the perturbation vector.
+func (p *SOMAT3APopulation) getPrtVector() []float64 {
+	prtv := make([]float64, p.Dimen)
+
+	for i := range prtv {
+		r := p.zeroToOneRNG.Rand()
+
+		switch {
+		case r < p.PRT:
+			prtv[i] = 1
+
+		case r > p.PRT:
+			prtv[i] = 0
+		}
+	}
+
+	return prtv
+}
+
+// getPrt calculates the adaptive PRT.
+func (p *SOMAT3APopulation) getPrt() float64 {
+	return 0.05 + 0.9*float64(p.fes/cec2020.GetMaxFES(p.Dimen))
+}
+
+// getStep calculates the step size.
+func (p *SOMAT3APopulation) getStep() float64 {
+	return 0.15 + 0.08*float64(p.fes/cec2020.GetMaxFES(p.Dimen))
+}

algo/ga/somat3aOrganisation.go

@@ -0,0 +1,76 @@
+// Copyright 2023 wanderer <a_mirre at utb dot cz>
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+package ga
+
+import "golang.org/x/exp/rand"
+
+func (p *SOMAT3APopulation) organise() (int, []int) {
+	leader := p.selectLeader()
+	migrants := p.selectMigrants(leader)
+
+	return leader, migrants
+}
+
+func (p *SOMAT3APopulation) selectMigrants(leader int) []int {
+	m := make([]int, p.M)
+	n := make([]int, p.N)
+	popCount := len(p.Population)
+
+	m[0] = rand.Intn(popCount)
+
+	for i := 1; i < p.M; i++ {
+	mgen:
+		candidateM := rand.Intn(popCount)
+
+		for _, v := range m {
+			if candidateM == v {
+				goto mgen
+			}
+		}
+
+		m[i] = candidateM
+	}
+
+	n[0] = rand.Intn(p.M)
+
+	for i := 1; i < p.N; i++ {
+	ngen:
+		candidateN := rand.Intn(p.M)
+
+		for _, v := range n {
+			if candidateN == v || candidateN == leader {
+				goto ngen
+			}
+		}
+
+		n[i] = candidateN
+	}
+
+	return n
+}
+
+func (p *SOMAT3APopulation) selectLeader() int {
+	k := make([]int, p.K)
+	popCount := len(p.Population)
+
+	k[0] = rand.Intn(popCount)
+
+	for i := 1; i < p.K; i++ {
+	kgen:
+		candidate := rand.Intn(popCount)
+
+		for _, v := range k {
+			if candidate == v {
+				goto kgen
+			}
+		}
+
+		k[i] = candidate
+	}
+
+	// choose one of K to become the leader.
+	leader := k[rand.Intn(p.K)]
+
+	return leader
+}

algo/ga/somat3aPopulation.go

@@ -0,0 +1,174 @@
+// Copyright 2023 wanderer <a_mirre at utb dot cz>
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+package ga
+
+import (
+	"git.dotya.ml/wanderer/math-optim/bench/cec2020"
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/stat/distuv"
+)
+
+// PopulationIndividual represents a single population individual.
+type SOMAT3APopulationIndividual struct {
+	CurX      DecisionVector
+	PRTVector []float64
+	Jumps     [][]float64
+}
+
+// ChampionIndividual is a representation of the best individual currently
+// available in the population.
+type SOMAT3AChampionIndividual struct {
+	X DecisionVector
+}
+
+// Population groups population individuals (agents) with metadata about the population.
+type SOMAT3APopulation struct {
+	// Population is a slice of population individuals.
+	Population []SOMAT3APopulationIndividual
+	// Champion represents the best individual of the population.
+	Champion SOMAT3AChampionIndividual
+	// Problem is the current benchmarking function this population is attempting to optimise.
+	Problem string
+	// ProblemFunction is the actual function to optimise.
+	ProblemFunc func([]float64) float64
+	// Dimen is the dimensionality of the problem being optimised.
+	Dimen int
+	// fes is the current number of function evaluations.
+	fes int
+	// Seed is the value used to (re)init population.
+	Seed         uint64
+	zeroToOneRNG distuv.Uniform
+
+	// PRT is the perturbation parameter.
+	PRT float64
+	// K is the number of individuals to choose the leader from.
+	K int
+	// M denotes how many individuals are picked from the population to form a
+	// `team` during the organisation phase, can be thought of as "team size".
+	M int
+	// N is the number of best individuals in each team selected for actual
+	// migration.
+	N int
+	// Njumps is the fixed number of jumps that each chosen migrating
+	// individual performs on their way to the leader.
+	Njumps int
+}
+
+// GetIndividual returns a reference to individual at position n.
+func (p *SOMAT3APopulation) GetIndividual(n uint) *PopulationIndividual {
+	return &PopulationIndividual{}
+}
+
+// GetBestIdx returns the index of the best population individual.
+func (p *SOMAT3APopulation) GetBestIdx() int {
+	f := p.ProblemFunc
+
+	bestIndividual := 0
+	// the first one is the best one.
+	bestVal := f(p.Population[0].CurX)
+
+	for i, v := range p.Population {
+		current := f(v.CurX)
+
+		if current < bestVal {
+			bestIndividual = i
+		}
+	}
+
+	return bestIndividual
+}
+
+// GetWorstIdx returns the index of the worst population individual.
+func (p *SOMAT3APopulation) GetWorstIdx() int {
+	f := p.ProblemFunc
+
+	worstIndividual := 0
+	// the first one is the worst one.
+	worstVal := f(p.Population[0].CurX)
+
+	for i, v := range p.Population {
+		current := f(v.CurX)
+
+		if current > worstVal {
+			worstIndividual = i
+		}
+	}
+
+	return worstIndividual
+}
+
+// Init initialises all individuals to random values.
+func (p *SOMAT3APopulation) Init() {
+	low := cec2020.SearchRange.Min()
+	high := cec2020.SearchRange.Max()
+	// zeroToOneRNG := distuv.Uniform{
+	p.zeroToOneRNG = distuv.Uniform{
+		Min: 0,
+		Max: 1,
+		Src: rand.NewSource(p.Seed),
+	}
+
+	// gaLogger.Printf("population initialisation - popCount: %d, seed: %d\n",
+	// len(p.Population), p.Seed,
+	// )
+
+	p.PRT = p.getPrt()
+
+	// for i,v  := range p.Population {
+	for i := range p.Population {
+		v := SOMAT3APopulationIndividual{}
+
+		v.CurX = make([]float64, p.Dimen)
+
+		for j := 0; j < p.Dimen; j++ {
+			// with current settings, this will always be `low + 0` since
+			// `high - low = 0`.
+			v.CurX[j] = low + p.zeroToOneRNG.Rand()*(high-low)
+		}
+
+		v.PRTVector = p.getPrtVector()
+
+		v.Jumps = make([][]float64, p.Njumps)
+
+		p.Population[i] = v
+	}
+
+	p.Champion = SOMAT3AChampionIndividual{X: p.Population[p.GetBestIdx()].CurX}
+}
+
+// Reinit reinitialises all individuals.
+func (p *SOMAT3APopulation) Reinit() {
+	p.Init()
+}
+
+// Clear sets all vectors to 0.
+func (p *SOMAT3APopulation) Clear() {
+	if p.Population != nil {
+		for _, v := range p.Population {
+			v.CurX = make([]float64, p.Dimen)
+		}
+	}
+}
+
+// Size returns the number of population individuals.
+func (p *SOMAT3APopulation) Size() int {
+	return len(p.Population)
+}
+
+// newPopulation returns a pointer to a new, uninitialised population.
+func newSOMAT3APopulation(benchProblem string, np, k, m, n, nj, dimen int) *SOMAT3APopulation {
+	p := &SOMAT3APopulation{
+		Problem:     benchProblem,
+		ProblemFunc: cec2020.Functions[benchProblem],
+		K:           k,
+		M:           m,
+		N:           n,
+		Njumps:      nj,
+		Dimen:       dimen,
+		Population:  make([]SOMAT3APopulationIndividual, np),
+	}
+
+	return p
+}