From 74109521b030c687f1f372aad8e65dc999450407 Mon Sep 17 00:00:00 2001 From: leo Date: Tue, 21 Feb 2023 23:51:52 +0100 Subject: [PATCH] ga: implement SOMA T3A algorithm --- algo/ga/error.go | 28 ++++ algo/ga/somat3a.go | 243 ++++++++++++++++++++++++++++++++ algo/ga/somat3aActualisation.go | 31 ++++ algo/ga/somat3aMigration.go | 80 +++++++++++ algo/ga/somat3aOrganisation.go | 76 ++++++++++ algo/ga/somat3aPopulation.go | 174 +++++++++++++++++++++++ 6 files changed, 632 insertions(+) create mode 100644 algo/ga/error.go create mode 100644 algo/ga/somat3aActualisation.go create mode 100644 algo/ga/somat3aMigration.go create mode 100644 algo/ga/somat3aOrganisation.go create mode 100644 algo/ga/somat3aPopulation.go diff --git a/algo/ga/error.go b/algo/ga/error.go new file mode 100644 index 0000000..ba7d2ef --- /dev/null +++ b/algo/ga/error.go @@ -0,0 +1,28 @@ +// Copyright 2023 wanderer +// 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} +} diff --git a/algo/ga/somat3a.go b/algo/ga/somat3a.go index 36b11db..30d891d 100644 --- a/algo/ga/somat3a.go +++ b/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{} +} diff --git a/algo/ga/somat3aActualisation.go b/algo/ga/somat3aActualisation.go new file mode 100644 index 0000000..9641d55 --- /dev/null +++ b/algo/ga/somat3aActualisation.go @@ -0,0 +1,31 @@ +// Copyright 2023 wanderer +// 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 + } +} diff --git a/algo/ga/somat3aMigration.go b/algo/ga/somat3aMigration.go new file mode 100644 index 0000000..a9c435e --- /dev/null +++ b/algo/ga/somat3aMigration.go @@ -0,0 +1,80 @@ +// Copyright 2023 wanderer +// 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)) +} diff --git a/algo/ga/somat3aOrganisation.go b/algo/ga/somat3aOrganisation.go new file mode 100644 index 0000000..b6f6816 --- /dev/null +++ b/algo/ga/somat3aOrganisation.go @@ -0,0 +1,76 @@ +// Copyright 2023 wanderer +// 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 +} diff --git a/algo/ga/somat3aPopulation.go b/algo/ga/somat3aPopulation.go new file mode 100644 index 0000000..6d523c2 --- /dev/null +++ b/algo/ga/somat3aPopulation.go @@ -0,0 +1,174 @@ +// Copyright 2023 wanderer +// 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 +}