diff --git a/algo/algo.go b/algo/algo.go index 4c296e5..39da143 100644 --- a/algo/algo.go +++ b/algo/algo.go @@ -9,20 +9,18 @@ import ( "sort" "sync" - "git.dotya.ml/wanderer/math-optim/algo/cec2020" "git.dotya.ml/wanderer/math-optim/algo/de" + "git.dotya.ml/wanderer/math-optim/algo/ga" "git.dotya.ml/wanderer/math-optim/bench" - c20 "git.dotya.ml/wanderer/math-optim/bench/cec2020" + "git.dotya.ml/wanderer/math-optim/bench/cec2020" "git.dotya.ml/wanderer/math-optim/report" "git.dotya.ml/wanderer/math-optim/stats" ) -// var Algos = []string{"Random Search", "Stochastic Hill Climbing"} - // mu protects access to meanStats. var mu sync.Mutex -// mCoMPL protexts access to comparisonOfMeansPicList. +// mCoMPL protects access to comparisonOfMeansPicList. var mCoMPL sync.Mutex var meanStats = &stats.MeanStats{} @@ -146,6 +144,95 @@ func PrepComparisonOfMeans(wg *sync.WaitGroup) (*report.PicList, int) { return pL, benchCount } +// PrepCEC2020ComparisonOfMeans prepares for comparison means of CEC2020 algos. +func PrepCEC2020ComparisonOfMeans(wg *sync.WaitGroup) (*report.PicList, int) { + pL := report.NewPicList() + meanStats := GetMeanStats() + algos := getAlgosFromAlgoMeans(meanStats.AlgoMeans) + + // construct title consisting of names of all involved algorithms. + for _, v := range algos { + switch pL.Algo { + case "": + pL.Algo = v + default: + pL.Algo += " vs " + v + } + } + + log.Println(`generating "Comparison of Means" plots`) + + algoCount := len(algos) + dimLen := len(cec2020.Dimensions) + benchCount := len(cec2020.Functions) + + for d := 0; d < dimLen; d++ { + // construct comparison for all benchmarking functions. + for i := 0; i < benchCount; i++ { + dimXAlgoMeanVals := make([]stats.AlgoMeanVals, 0, algoCount) + + for j := 0; j < algoCount; j++ { + ms := &stats.AlgoMeanVals{ + Title: meanStats.AlgoMeans[i+(j*benchCount)].Algo, + MeanVals: meanStats.AlgoMeans[i+(j*benchCount)].BenchMeans[d].MeanVals, + } + + dimXAlgoMeanVals = append(dimXAlgoMeanVals, *ms) + } + + dimens := meanStats.AlgoMeans[i].BenchMeans[d].Dimens + iterations := meanStats.AlgoMeans[i].BenchMeans[d].Iterations + bench := meanStats.AlgoMeans[i].BenchMeans[d].Bench + + wg.Add(1) + + // construct plots concurrently. + go PlotMeanValsMulti( + wg, dimens, iterations, bench, "plot-", ".pdf", + dimXAlgoMeanVals..., + ) + } + } + + // wait for all plotting goroutines. + wg.Wait() + + pL.Pics = getComparisonOfMeansPics() + + return pL, dimLen +} + +// getAlgosFromAlgoMeans extracts algorithms used from the means list and +// returns it as a []string. +func getAlgosFromAlgoMeans(s []stats.AlgoBenchMean) []string { + algos := make([]string, 0) + + // learn how many algos were processed based on the data. + for _, v := range s { + // if algos is empty just add the value directly, else determine if + // it's already been added or not. + if len(algos) > 0 { + alreadyadded := false + + for _, algoName := range algos { + if algoName == v.Algo { + // early bail if already added. + alreadyadded = true + break + } + } + + if !alreadyadded { + algos = append(algos, v.Algo) + } + } else { + algos = append(algos, v.Algo) + } + } + + return algos +} + // DoRandomSearch executes a search using the 'Random search' method. func DoRandomSearch(wg *sync.WaitGroup, m *sync.Mutex) { defer wg.Done() @@ -418,10 +505,8 @@ func DojDE(wg *sync.WaitGroup, m *sync.Mutex) { func DoCEC2020jDE(wg *sync.WaitGroup, m *sync.Mutex) { defer wg.Done() - cec2020.LogPrintln("starting") - // funcCount is the number of bench functions available and tested. - funcCount := len(c20.Functions) + funcCount := len(cec2020.Functions) // stats for the current algo. algoStats := make([][]stats.Stats, funcCount) // ch serves as a way to get the actual computed output. @@ -438,7 +523,7 @@ func DoCEC2020jDE(wg *sync.WaitGroup, m *sync.Mutex) { cr := 0.9 for i := range algoStats { - jDE := cec2020.NewjDE() + jDE := ga.NewjDE() // params: // Generations, minimum bench iterations, mutation strategy, parameter @@ -451,7 +536,7 @@ func DoCEC2020jDE(wg *sync.WaitGroup, m *sync.Mutex) { // 0..17 to choose a mutation strategy, // 0..1 to select a parameter self-adaptation scheme, // np >= 4 as initial population size. - jDE.Init(-1, 30, 0, 0, np, f, cr, c20.Dimensions, c20.FuncNames[i], ch, chAlgoMeans) + jDE.Init(-1, 30, 0, 0, np, f, cr, cec2020.Dimensions, cec2020.FuncNames[i], ch, chAlgoMeans) go jDE.Run() } @@ -466,8 +551,8 @@ func DoCEC2020jDE(wg *sync.WaitGroup, m *sync.Mutex) { saveAlgoMeans(*aM) } - pCh := make(chan report.PicList, funcCount*len(c20.Dimensions)) - pMeanCh := make(chan report.PicList, funcCount*len(c20.Dimensions)) + pCh := make(chan report.PicList, funcCount*len(cec2020.Dimensions)) + pMeanCh := make(chan report.PicList, funcCount*len(cec2020.Dimensions)) for _, algoStat := range algoStats { go plotAllDims(algoStat, "plot", ".pdf", pCh, pMeanCh) @@ -493,3 +578,87 @@ func DoCEC2020jDE(wg *sync.WaitGroup, m *sync.Mutex) { stats.SaveTable(algoName, algoStats) m.Unlock() } + +// DoCEC2020SOMAT3A performs a search using the SOMA T3A method. +func DoCEC2020SOMAT3A(wg *sync.WaitGroup, m *sync.Mutex) { + defer wg.Done() + + // funcCount is the number of bench functions available and tested. + funcCount := len(cec2020.Functions) + // stats for the current algo. + algoStats := make([][]stats.Stats, funcCount) + // ch serves as a way to get the actual computed output. + ch := make(chan []stats.Stats, funcCount) + chAlgoMeans := make(chan *stats.AlgoBenchMean, funcCount) + + defer close(ch) + defer close(chAlgoMeans) + + // somat3a params. + np := 50 + k := 10 + mSize := 10 + n := 5 + njumps := 7 + + for i := range algoStats { + somat3a := ga.NewSOMAT3A() + + // params: + // Generations, minimum bench iterations, initial population size, + // leader candidates, migration candidates group size, number of + // migrants, number of jumps each migrant performs, dimensions, bench + // name and synchronisation channels. + // + // -1 to disable generation limits, + // n > 0 for minimum bench iterations, + // np >= k+mSize as initial population size, + // k, mSize, n and njumps >= 0. + err := somat3a.Init(-1, 30, np, k, mSize, n, njumps, + cec2020.Dimensions, cec2020.FuncNames[i], + ch, chAlgoMeans, + ) + if err != nil { + log.Panicf("Failed to initialise SOMA T3A, error: %q", err) + } + + go somat3a.Run() + } + + // get results. + for i := range algoStats { + s := <-ch + aM := <-chAlgoMeans + + algoStats[i] = s + + saveAlgoMeans(*aM) + } + + pCh := make(chan report.PicList, funcCount*len(cec2020.Dimensions)) + pMeanCh := make(chan report.PicList, funcCount*len(cec2020.Dimensions)) + + for _, algoStat := range algoStats { + go plotAllDims(algoStat, "plot", ".pdf", pCh, pMeanCh) + } + + pLs := []report.PicList{} + pLsMean := []report.PicList{} + + for range algoStats { + pL := <-pCh + pLMean := <-pMeanCh + + pLs = append(pLs, pL) + pLsMean = append(pLsMean, pLMean) + } + + algoName := "SOMA T3A" + + // protect access to shared data. + m.Lock() + report.SavePicsToFile(pLs, pLsMean, algoName) + + stats.SaveTable(algoName, algoStats) + m.Unlock() +} diff --git a/algo/cec2020/doc.go b/algo/cec2020/doc.go deleted file mode 100644 index 3a174eb..0000000 --- a/algo/cec2020/doc.go +++ /dev/null @@ -1,5 +0,0 @@ -// Copyright 2023 wanderer -// SPDX-License-Identifier: GPL-3.0-or-later - -// Package cec2020 implements algorithms required by CEC2020: jDE and T3A. -package cec2020 diff --git a/algo/cec2020/jDE.go b/algo/cec2020/jDE.go deleted file mode 100644 index 9ab4c73..0000000 --- a/algo/cec2020/jDE.go +++ /dev/null @@ -1,500 +0,0 @@ -// Copyright 2023 wanderer -// SPDX-License-Identifier: GPL-3.0-or-later - -package cec2020 - -import ( - "sort" - "time" - - "golang.org/x/exp/rand" - - "git.dotya.ml/wanderer/math-optim/bench/cec2020" - "git.dotya.ml/wanderer/math-optim/stats" - "gonum.org/v1/gonum/stat/distuv" -) - -// JDE is a holder for the settings of an instance of a self-adapting -// differential evolution (jDE) algorithm. -type JDE struct { - // Generations denotes the number of generations the population evolves - // for. Special value -1 disables limiting the number of generations. - Generations int - // BenchMinIters is the number of iterations the bench function will be re-run (for statistical purposes). - BenchMinIters int - // Dimensions to solve the problem for. - Dimensions []int - // F is the initial value of the differential weight (mutation/weighting factor). - F float64 - // CR is the initial value of the crossover probability constant. - CR float64 - // rngF is a random number generator for differential weight (mutation/weighting factor). - rngF distuv.Uniform - // cr is a random number generator for crossover probability constant adapted over time. - rngCR distuv.Uniform - // MutationStrategy selects the mutation strategy, i.e. the variant of the - // jDE algorithm (0..17), see mutationStrategies.go for more details. - MutationStrategy int - // AdptScheme is the parameter self-adaptation scheme (0..1). - AdptScheme int - // NP is the initial population size. - NP int - // BenchName is a name of the problem to optimise. - BenchName string - // ch is a channel for writing back computed results. - ch chan []stats.Stats - // chAlgoMeans is a channel for writing back algo means. - chAlgoMeans chan *stats.AlgoBenchMean - // initialised denotes the initialisation state of the struct. - initialised bool -} - -const ( - // jDEMinNP is the minimum size of the initial population for jDE. - // for jDE PaGMO specifies 8. - jDEMinNP = 4 - // fMin is the minimum allowed value of the differential weight. - fMin = 0.36 - // fMax is the maximum allowed value of the differential weight. - fMax = 1.0 - // crMin is the minimum allowed value of the crossover probability constant. - crMin = 0.0 - // crMax is the maximum allowed value of the crossover probability constant. - crMax = 1.0 - - // tau1 is used in parameter self-adaptation. - tau1 = 0.1 - // tau2 is used in parameter self-adaptation. - tau2 = 0.1 - // fl is used in parameter self-adaptation. - fl = 0.1 - // fu is used in parameter self-adaptation. - fu = 0.9 -) - -// Init initialises the jDE algorithm, performs sanity checks on the inputs. -func (j *JDE) Init(generations, benchMinIters, mutStrategy, adptScheme, np int, f, cr float64, dimensions []int, bench string, ch chan []stats.Stats, chAlgoMeans chan *stats.AlgoBenchMean) { - if j == nil { - cec2020Logger.Fatalln("jDE needs to be initialised before calling RunjDE, exiting...") - } - - // check input parameters. - switch { - case generations == 0: - cec2020Logger.Fatalln("Generations cannot be 0, got", generations) - - case generations == -1: - cec2020Logger.Println("Generations is '-1', disabling generation limits..") - - case benchMinIters < 1: - cec2020Logger.Fatalln("Minimum bench iterations cannot be less than 1, got:", benchMinIters) - - case mutStrategy < 0 || mutStrategy > 17: - cec2020Logger.Fatalln("Mutation strategy needs to be from the interval <0; 17>, got", mutStrategy) - - case adptScheme < 0 || adptScheme > 1: - cec2020Logger.Fatalln("Parameter self-adaptation scheme needs to be from the interval <0; 1>, got", adptScheme) - - case np < jDEMinNP: - cec2020Logger.Fatalf("NP cannot be less than %d, got: %d\n.", jDEMinNP, np) - - case f < fMin || f > fMax: - cec2020Logger.Fatalf("F needs to be from the interval <%f;%f>, got: %f\n.", fMin, fMax, f) - - case cr < crMin || cr > crMax: - cec2020Logger.Fatalf("CR needs to be from the interval <%f;>%f, got: %f\n.", crMin, crMax, cr) - - case len(dimensions) == 0: - cec2020Logger.Fatalf("Dimensions cannot be empty, got: %+v\n", dimensions) - - case bench == "": - cec2020Logger.Fatalln("Bench cannot be empty, got:", bench) - } - - j.Generations = generations - j.BenchMinIters = benchMinIters - j.MutationStrategy = mutStrategy - j.AdptScheme = adptScheme - j.NP = np - j.F = f - j.CR = cr - - rngsrc := rand.NewSource(uint64(rand.Int63())) - - j.rngF = distuv.Uniform{Min: fMin, Max: fMax, Src: rngsrc} - j.rngCR = distuv.Uniform{Min: crMin, Max: crMax, Src: rngsrc} - - j.Dimensions = dimensions - j.BenchName = bench - j.ch = ch - j.chAlgoMeans = chAlgoMeans - - j.initialised = true - - cec2020Logger.Printf("jDE init done, jDE:%+v", j) -} - -// InitAndRun initialises the jDE algorithm, performs sanity checks on the -// inputs and calls the Run method. -func (j *JDE) InitAndRun(generations, benchMinIters, mutStrategy, adptScheme, np int, f, cr float64, dimensions []int, bench string, ch chan []stats.Stats, chAlgoMeans chan *stats.AlgoBenchMean) { - if j == nil { - cec2020Logger.Fatalln("jDE is nil, NewjDE() needs to be called first. exiting...") - } - - j.Init(generations, benchMinIters, mutStrategy, adptScheme, np, f, cr, dimensions, bench, ch, chAlgoMeans) - - j.Run() -} - -// Run self-adapting differential evolution algorithm. -func (j *JDE) Run() { - if j == nil { - cec2020Logger.Fatalln("jDE is nil, NewjDE() needs to be called first. exiting...") - } - - if !j.initialised { - cec2020Logger.Fatalln("jDE needs to be initialised before calling Run(), exiting...") - } - - var jDEStats []stats.Stats - - jDEMeans := &stats.AlgoBenchMean{ - Algo: "jDE", - BenchMeans: make([]stats.BenchMean, 0, len(j.Dimensions)), - } - - // run evolve for all dimensions. - for _, dim := range j.Dimensions { - maxFES := cec2020.GetMaxFES(dim) - fesPerIter := int(float64(maxFES / j.NP)) - jDEStatDimX := &stats.Stats{ - Algo: "jDE", - Dimens: dim, - Iterations: j.BenchMinIters, - Generations: maxFES, - NP: j.NP, - F: j.F, - CR: j.CR, - } - funcStats := &stats.FuncStats{BenchName: j.BenchName} - dimXMean := &stats.BenchMean{ - Bench: j.BenchName, - Dimens: dim, - Iterations: j.BenchMinIters, - Generations: maxFES, - Neighbours: -1, - } - uniDist := distuv.Uniform{ - Min: cec2020.SearchRange.Min(), - Max: cec2020.SearchRange.Max(), - } - - if maxFES == -1 { - cec2020Logger.Fatalf("could not get maxFES for current dim (%d), bailing", dim) - } - - cec2020Logger.Printf("running bench \"%s\" for %dD, maxFES: %d\n", - j.BenchName, dim, maxFES, - ) - - funcStats.BenchResults = make([]stats.BenchRound, j.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 < j.BenchMinIters; iter++ { - cec2020Logger.Printf("run: %d, bench: %s, %dD, started at %s", - iter, j.BenchName, dim, start, - ) - - funcStats.BenchResults[iter].Iteration = iter - - uniDist.Src = src - - // create a population with known params. - pop := newPopulation(j.BenchName, j.NP, dim) - - // set population seed. - pop.Seed = uint64(time.Now().UnixNano()) - // initialise the population. - pop.Init() - - var bestResult float64 - - // the core. - for i := 0; i < fesPerIter; i++ { - r := j.evaluate(pop) - - // distinguish the first or any of the subsequent iterations. - switch i { - case 0: - bestResult = r - - default: - // call evolve where jDE runs. - j.evolve(pop, &uniDist) - - // take measurements of current population fitness. - r = j.evaluate(pop) - - // save if better. - if r < bestResult { - bestResult = r - } - } - - funcStats.BenchResults[iter].Results = append( - funcStats.BenchResults[iter].Results, - bestResult, - ) - // 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 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 < j.NP-1; x++ { - funcStats.BenchResults[iter].Results = append( - funcStats.BenchResults[iter].Results, - bestResult, - ) - } - } - } - - elapsed := time.Since(start) - - cec2020Logger.Printf("completed: bench: %s, %dD, computing took %s\n", - j.BenchName, dim, elapsed, - ) - - // get mean vals. - dimXMean.MeanVals = stats.GetMeanVals(funcStats.BenchResults, maxFES) - funcStats.MeanVals = dimXMean.MeanVals - - jDEStatDimX.BenchFuncStats = append( - jDEStatDimX.BenchFuncStats, *funcStats, - ) - - jDEStats = append(jDEStats, *jDEStatDimX) - - jDEMeans.BenchMeans = append(jDEMeans.BenchMeans, *dimXMean) - } - - sort.Sort(jDEMeans) - - j.chAlgoMeans <- jDEMeans - j.ch <- jDEStats -} - -// evaluate evaluates the fitness of current population. -func (j *JDE) evaluate(pop *Population) float64 { - f := cec2020.Functions[pop.Problem] - - bestIndividual := pop.Population[pop.GetBestIdx()].CurX - bestSolution := f(bestIndividual) - - return bestSolution -} - -// evolve evolves a population by running the jDE (self-adapting Differential -// Evolution) algorithm on the passed population until termination conditions -// are met. -// nolint: gocognit -func (j *JDE) evolve(pop *Population, uniDist *distuv.Uniform) { - popCount := len(pop.Population) - - for i, currentIndividual := range pop.Population { - idcs := make([]int, 0, popCount-1) - // idcs := make([]int, popCount-1) - - // gather indices. - for k := 0; k < popCount; k++ { - if k != i { - idcs = append(idcs, k) - } - } - - // randomly choose 3 of those idcs. - selectedIdcs := make([]int, 0) - - selectedA := false - selectedB := false - selectedC := false - - for !selectedA { - candidateA := rand.Intn(len(idcs)) % len(idcs) - - if candidateA != i { - selectedIdcs = append(selectedIdcs, candidateA) - selectedA = true - } - } - - for !selectedB { - a := selectedIdcs[0] - candidateB := rand.Intn(len(idcs)) % len(idcs) - - if candidateB != i && candidateB != a { - selectedIdcs = append(selectedIdcs, candidateB) - selectedB = true - } - } - - for !selectedC { - a := selectedIdcs[0] - b := selectedIdcs[1] - candidateC := rand.Intn(len(idcs)) % len(idcs) - - if candidateC != i && candidateC != a && candidateC != b { - selectedIdcs = append(selectedIdcs, candidateC) - selectedC = true - } - } - - // selected contains the selected population individuals. - selected := make([]PopulationIndividual, 0) - - // select individuals for rand/1/bin - for _, idx := range selectedIdcs { - for k := 0; k < popCount; k++ { - if k == idx { - selected = append(selected, pop.Population[idx]) - } - } - } - - mutant := make([]float64, pop.Dimen) - - // mutate. - for k := 0; k < pop.Dimen; k++ { - // mutant = a + mut * (b - c) - mutant[k] = selected[0].CurX[k] + (pop.BestF * (selected[1].CurX[k] - selected[2].CurX[k])) - - // clip values to <0;1>. - switch { - case mutant[k] < 0: - mutant[k] = 0 - - case mutant[k] > 1: - mutant[k] = 1 - } - } - - crossPoints := make([]bool, pop.Dimen) - - // prepare crossover points (binomial crossover). - for k := 0; k < pop.Dimen; k++ { - if v := uniDist.Rand(); v < pop.BestCR { - crossPoints[k] = true - } else { - crossPoints[k] = false - } - } - - trial := make([]float64, pop.Dimen) - - // recombine using crossover points. - for k := 0; k < pop.Dimen; k++ { - if crossPoints[k] { - trial[k] = currentIndividual.CurX[k] - } else { - trial[k] = mutant[k] - } - } - - f := cec2020.Functions[pop.Problem] - currentFitness := f(currentIndividual.CurX) - trialFitness := f(trial) - - // replace if better. - if trialFitness < currentFitness { - pop.Population[i].CurX = trial - } - } - - // adapt parameters. - j.adaptParameters(pop) -} - -// adptParameters adapts parameters F and CR based on -// https://labraj.feri.um.si/images/0/05/CEC09_slides_Brest.pdf. -func (j *JDE) adaptParameters(p *Population) { - var nuF float64 - - var nuCR float64 - - switch { - case j.AdptScheme == 0: - if rand2 := j.getRandF(); rand2 < tau1 { - nuF = fl + (j.getRandF() * fu) - } else { - nuF = p.f[len(p.f)-1] - } - - if rand4 := j.getRandCR(); rand4 < tau2 { - nuCR = j.getRandCR() - } else { - nuCR = p.cr[len(p.cr)-1] - } - - case j.AdptScheme == 1: - nuF = p.BestF + j.getRandF()*0.5 - nuCR = p.BestCR + j.getRandCR()*0.5 - } - - // sort out F. - p.f = append(p.f, nuF) - // sort out CR - p.cr = append(p.cr, nuCR) - - // update best, if improved. - switch { - case nuF < p.BestF: - p.BestF = nuF - - case nuCR < p.BestCR: - p.BestCR = nuCR - } -} - -// getRandF returns a random value of F for parameter self-adaptation. -func (j *JDE) getRandF() float64 { - return j.rngF.Rand() -} - -// getRandCR returns a random value of CR for parameter self-adaptation. -func (j *JDE) getRandCR() float64 { - return j.rngCR.Rand() -} - -// NewjDE returns a pointer to a new, uninitialised jDE instance. -func NewjDE() *JDE { - return &JDE{} -} - -// LogPrintln wraps the jDE logger's Println func. -func LogPrintln(v ...any) { - cec2020Logger.Println(v...) -} - -// LogPrintf wraps the jDE logger's Printf func. -func LogPrintf(s string, v ...any) { - cec2020Logger.Printf(s, v...) -} - -// LogFatalln wraps the jDE logger's Fatalln func. -func LogFatalln(s string) { - cec2020Logger.Fatalln(s) -} - -// LogFatalf wraps the jDE logger's Fatalf func. -func LogFatalf(s string, v ...any) { - cec2020Logger.Fatalf(s, v...) -} diff --git a/algo/cec2020/log.go b/algo/cec2020/log.go deleted file mode 100644 index d6d1cb7..0000000 --- a/algo/cec2020/log.go +++ /dev/null @@ -1,12 +0,0 @@ -// Copyright 2023 wanderer -// SPDX-License-Identifier: GPL-3.0-or-later - -package cec2020 - -import ( - "log" - "os" -) - -// cec2020Logger declares and initialises a "custom" CEC2020 logger. -var cec2020Logger = log.New(os.Stderr, " *** ∁ cec2020:", log.Ldate|log.Ltime|log.Lshortfile) diff --git a/algo/cec2020/population.go b/algo/cec2020/population.go deleted file mode 100644 index 0af1931..0000000 --- a/algo/cec2020/population.go +++ /dev/null @@ -1,167 +0,0 @@ -// Copyright 2023 wanderer -// SPDX-License-Identifier: GPL-3.0-or-later - -package cec2020 - -import ( - "git.dotya.ml/wanderer/math-optim/bench/cec2020" - "golang.org/x/exp/rand" - - "gonum.org/v1/gonum/stat/distuv" -) - -type ( - // DecisionVector is a []float64 abstraction representing the decision vector. - DecisionVector []float64 - // FitnessVector is a []float64 abstraction representing the fitness vector. - FitnessVector []float64 - // ConstraintVector is a []float64 abstraction representing the constraint vector. - ConstraintVector []float64 -) - -// PopulationIndividual representats a single population individual. -type PopulationIndividual struct { - CurX DecisionVector - CurV DecisionVector - CurC ConstraintVector - CurF FitnessVector - - BestX DecisionVector - BestC ConstraintVector - BestF FitnessVector -} - -// ChampionIndividual is a representation of the best individual currently -// available in the population. -type ChampionIndividual struct { - X DecisionVector - C ConstraintVector - F FitnessVector -} - -// Population groups population individuals (agents) with metadata about the population. -type Population struct { - // Population is a slice of population individuals. - Population []PopulationIndividual - // Problem is the current benchmarking function this population is attempting to optimise. - Problem string - // Dimen is the dimensionality of the problem being optimised. - Dimen int - // Seed is the value used to (re)init population. - Seed uint64 - - // f is the differential weight (mutation/weighting factor) adapted over time. - f []float64 - // cr is the crossover probability constant adapted over time. - cr []float64 - // BestF is the best recorded value of the differential weight F. - BestF float64 - // BestCR is the best recorded value of the differential weight CR. - BestCR float64 -} - -// GetIndividual returns a reference to individual at position n. -func (p *Population) GetIndividual(n uint) *PopulationIndividual { return &PopulationIndividual{} } - -// GetBestIdx returns the index of the best population individual. -func (p *Population) GetBestIdx() int { - f := cec2020.Functions[p.Problem] - - 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 *Population) GetWorstIdx() int { - f := cec2020.Functions[p.Problem] - - 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 *Population) Init() { - uniform := distuv.Uniform{ - Min: cec2020.SearchRange.Min(), - Max: cec2020.SearchRange.Max(), - } - uniform.Src = rand.NewSource(p.Seed) - - for i, v := range p.Population { - v.CurX = make([]float64, p.Dimen) - - for j := 0; j < p.Dimen; j++ { - v.CurX[j] = uniform.Rand() - } - - p.Population[i] = v - } - - p.f = make([]float64, p.Size()) - p.cr = make([]float64, p.Size()) -} - -// Reinit reinitialises all individuals. -func (p *Population) Reinit() { - p.Init() -} - -// ReinitN reinitialises the individual at position n. -func (p *Population) ReinitN(n uint) {} - -// Clear sets all vectors to 0. -func (p *Population) Clear() { - if p.Population != nil { - for _, v := range p.Population { - v.CurX = make([]float64, p.Dimen) - v.CurC = make([]float64, p.Dimen) - v.CurF = make([]float64, p.Dimen) - v.CurV = make([]float64, p.Dimen) - v.BestX = make([]float64, p.Dimen) - v.BestC = make([]float64, p.Dimen) - v.BestF = make([]float64, p.Dimen) - } - } -} - -// MeanVelocity computes the mean current velocity of all individuals in the population. -func (p *Population) MeanVelocity() float64 { return 0.0 } - -// Size returns the number of population individuals. -func (p *Population) Size() int { return len(p.Population) } - -// newPopulation returns a pointer to a new, uninitialised population. -func newPopulation(benchProblem string, np, dimen int) *Population { - p := &Population{} - - p.Problem = benchProblem - - p.Dimen = dimen - - // pre-alloc. - p.Population = make([]PopulationIndividual, np) - - return p -} diff --git a/algo/ga/de.go b/algo/ga/de.go new file mode 100644 index 0000000..d74d4b7 --- /dev/null +++ b/algo/ga/de.go @@ -0,0 +1,411 @@ +// Copyright 2023 wanderer +// SPDX-License-Identifier: GPL-3.0-or-later + +package ga + +import ( + "sort" + "time" + + "golang.org/x/exp/rand" + + "git.dotya.ml/wanderer/math-optim/bench" + "git.dotya.ml/wanderer/math-optim/stats" + "gonum.org/v1/gonum/stat/distuv" +) + +// DE is a holder for the settings of an instance of a Differential Evolution +// (DE) algorithm. +type DE struct { + // Generations denotes the number of generations the population evolves + // for. Special value -1 disables limiting the number of generations. + Generations int + // BenchMinIters is the number of iterations the bench function will be re-run (for statistical purposes). + BenchMinIters int + // Dimensions to solve the problem for. + Dimensions []int + // F is the differential weight (mutation/weighting factor). + F float64 + // CR is the crossover probability constant. + CR float64 + // MutationStrategy selects the mutation strategy, i.e. the variant of the + // DE algorithm (0..17), see mutationStrategies.go for more details. + MutationStrategy int + // NP is the initial population size. + NP int + // BenchName is a name of the problem to optimise. + BenchName string + // ch is a channel for writing back computed results. + ch chan []stats.Stats + // chAlgoMeans is a channel for writing back algo means. + chAlgoMeans chan *stats.AlgoBenchMean + // initialised denotes the initialisation state of the struct. + initialised bool +} + +const ( + // MinNPDE is the minimum size of the initial population for DE. + minNPDE = 4 + // fMin is the minimum allowed value of the differential weight. + fMinDE = 0.5 + // fMax is the maximum allowed value of the differential weight. + fMaxDE = 2.0 + // crMin is the minimum allowed value of the crossover probability constant. + crMinDE = 0.2 + // crMax is the maximum allowed value of the crossover probability constant. + crMaxDE = 0.9 +) + +// dELogger is a "custom" DE logger. +var dELogger = newLogger(" *** δ DE:") + +// Init initialises the DE algorithm, performs sanity checks on the inputs. +func (d *DE) Init( + generations, benchMinIters, mutStrategy, np int, + f, cr float64, + dimensions []int, + bench string, + ch chan []stats.Stats, + chAlgoMeans chan *stats.AlgoBenchMean, +) { + if d == nil { + dELogger.Fatalln("DE needs to be initialised before calling Run(), exiting...") + } + + // check input parameters. + switch { + case generations == 0: + dELogger.Fatalln("Generations cannot be 0, got", generations) + + case generations == -1: + dELogger.Println("Generations is '-1', disabling generation limits..") + + case benchMinIters < 1: + dELogger.Fatalln("Minimum bench iterations cannot be less than 1, got:", benchMinIters) + + case mutStrategy < 0 || mutStrategy > 17: + dELogger.Fatalln("Mutation strategy needs to be from the interval <0; 17>, got", mutStrategy) + + case np < minNPDE: + dELogger.Fatalf("NP cannot be less than %d, got: %d\n.", minNPDE, np) + + case f < fMinDE || f > fMaxDE: + dELogger.Fatalf("F needs to be from the interval <%f;%f>, got: %f\n.", fMinDE, fMaxDE, f) + + case cr < crMinDE || cr > crMaxDE: + dELogger.Fatalf("CR needs to be from the interval <%f;>%f, got: %f\n.", crMinDE, crMaxDE, cr) + + case len(dimensions) == 0: + dELogger.Fatalf("Dimensions cannot be empty, got: %+v\n", dimensions) + + case bench == "": + dELogger.Fatalln("Bench cannot be unset, got:", bench) + } + + d.Generations = generations + d.BenchMinIters = benchMinIters + d.MutationStrategy = mutStrategy + d.NP = np + d.F = f + d.CR = cr + d.Dimensions = dimensions + d.BenchName = bench + d.ch = ch + d.chAlgoMeans = chAlgoMeans + + d.initialised = true +} + +// InitAndRun initialises the DE algorithm, performs sanity checks on the +// inputs and calls the Run method. +func (d *DE) InitAndRun( + generations, benchMinIters, mutStrategy, np int, + f, cr float64, + dimensions []int, + bench string, + ch chan []stats.Stats, + chAlgoMeans chan *stats.AlgoBenchMean, +) { + if d == nil { + dELogger.Fatalln("DE is nil, NewDE() needs to be called first. exiting...") + } + + d.Init( + generations, benchMinIters, mutStrategy, np, + f, cr, + dimensions, + bench, + ch, + chAlgoMeans, + ) + + d.Run() +} + +// Run self-adapting differential evolution algorithm. +func (d *DE) Run() { + if d == nil { + dELogger.Fatalln("DE is nil, NewDE() needs to be called first. exiting...") + } + + if !d.initialised { + dELogger.Fatalln("DE needs to be initialised before calling Run(), exiting...") + } + + var dEStats []stats.Stats + + dEMeans := &stats.AlgoBenchMean{ + Algo: "DE", + BenchMeans: make([]stats.BenchMean, 0, len(d.Dimensions)), + } + + // run evolve for for all dimensions. + for _, dim := range d.Dimensions { + maxFES := bench.GetGAMaxFES(dim) + fesPerIter := int(float64(maxFES / d.NP)) + dEStatDimX := &stats.Stats{ + Algo: "DE", + Dimens: dim, + Iterations: d.BenchMinIters, + Generations: maxFES, + NP: d.NP, + F: d.F, + CR: d.CR, + } + funcStats := &stats.FuncStats{BenchName: d.BenchName} + dimXMean := &stats.BenchMean{ + Bench: d.BenchName, + Dimens: dim, + Iterations: d.BenchMinIters, + Generations: maxFES, + Neighbours: -1, + } + benchFuncParams := bench.FunctionParams[d.BenchName] + uniDist := distuv.Uniform{ + Min: benchFuncParams.Min(), + Max: benchFuncParams.Max(), + } + + funcStats.BenchResults = make([]stats.BenchRound, d.BenchMinIters) + + // rand.Seed(uint64(time.Now().UnixNano())) + // create and seed a source of preudo-randomness + src := rand.NewSource(uint64(rand.Int63())) + + // track execution duration. + start := time.Now() + + for iter := 0; iter < d.BenchMinIters; iter++ { + dELogger.Printf("run: %d, bench: %s, %dD, started at %s", + iter, d.BenchName, dim, start, + ) + + funcStats.BenchResults[iter].Iteration = iter + + uniDist.Src = src + + // create a population with known params. + pop := newPopulation(d.BenchName, d.NP, dim) + + // set population seed. + pop.Seed = uint64(time.Now().UnixNano()) + // initialise the population. + pop.Init() + + var bestResult float64 + + // the core. + for i := 0; i < fesPerIter; i++ { + r := d.evaluate(pop) + + // distinguish the first or any of the subsequent iterations. + switch i { + case 0: + bestResult = r + + default: + // call evolve where actual DE runs. + d.evolve(pop, &uniDist) + + // take measurements of current population fitness. + r = d.evaluate(pop) + + // save if better. + if r < bestResult { + bestResult = r + } + } + + funcStats.BenchResults[iter].Results = append( + funcStats.BenchResults[iter].Results, + bestResult, + ) + // 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 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 < d.NP-1; x++ { + funcStats.BenchResults[iter].Results = append( + funcStats.BenchResults[iter].Results, + bestResult, + ) + } + } + } + + elapsed := time.Since(start) + + dELogger.Printf("completed: bench: %s, %dD, computing took %s\n", + d.BenchName, dim, elapsed, + ) + + // get mean vals. + dimXMean.MeanVals = stats.GetMeanVals(funcStats.BenchResults, maxFES) + funcStats.MeanVals = dimXMean.MeanVals + + dEStatDimX.BenchFuncStats = append( + dEStatDimX.BenchFuncStats, *funcStats, + ) + + dEStats = append(dEStats, *dEStatDimX) + + dEMeans.BenchMeans = append(dEMeans.BenchMeans, *dimXMean) + } + + sort.Sort(dEMeans) + + d.chAlgoMeans <- dEMeans + d.ch <- dEStats +} + +// evaluate evaluates the fitness of current population. +func (d *DE) evaluate(pop *Population) float64 { + f := bench.Functions[pop.Problem] + + bestIndividual := pop.Population[pop.GetBestIdx()].CurX + bestSolution := f(bestIndividual) + + return bestSolution +} + +// evolve evolves a population by running the DE (Differential Evolution) +// algorithm on the passed population until termination conditions are met. +// nolint: gocognit +func (d *DE) evolve(pop *Population, uniDist *distuv.Uniform) { + popCount := len(pop.Population) + + for i, currentIndividual := range pop.Population { + idcs := make([]int, 0, popCount-1) + + // gather indices. + for k := 0; k < popCount; k++ { + if k != i { + idcs = append(idcs, k) + } + } + + // randomly choose 3 of those idcs. + selectedIdcs := make([]int, 0) + + selectedA := false + selectedB := false + selectedC := false + + for !selectedA { + candidateA := rand.Intn(len(idcs)) % len(idcs) + + if candidateA != i { + selectedIdcs = append(selectedIdcs, candidateA) + selectedA = true + } + } + + for !selectedB { + a := selectedIdcs[0] + candidateB := rand.Intn(len(idcs)) % len(idcs) + + if candidateB != i && candidateB != a { + selectedIdcs = append(selectedIdcs, candidateB) + selectedB = true + } + } + + for !selectedC { + a := selectedIdcs[0] + b := selectedIdcs[1] + candidateC := rand.Intn(len(idcs)) % len(idcs) + + if candidateC != i && candidateC != a && candidateC != b { + selectedIdcs = append(selectedIdcs, candidateC) + selectedC = true + } + } + + // selected contains the selected population individuals. + selected := make([]PopulationIndividual, 0) + + // select individuals for rand/1/bin + for _, idx := range selectedIdcs { + for k := 0; k < popCount; k++ { + if k == idx { + selected = append(selected, pop.Population[idx]) + } + } + } + + mutant := make([]float64, pop.Dimen) + + // mutate. + for k := 0; k < pop.Dimen; k++ { + // mutant = a + mut * (b - c) + mutant[k] = selected[0].CurX[k] + (d.F * (selected[1].CurX[k] - selected[2].CurX[k])) + + // clip values to <0;1>. + switch { + case mutant[k] < 0: + mutant[k] = 0 + + case mutant[k] > 1: + mutant[k] = 1 + } + } + + crossPoints := make([]bool, pop.Dimen) + + // prepare crossover points (binomial crossover). + for k := 0; k < pop.Dimen; k++ { + if v := uniDist.Rand(); v < d.CR { + crossPoints[k] = true + } else { + crossPoints[k] = false + } + } + + trial := make([]float64, pop.Dimen) + + // recombine using crossover points. + for k := 0; k < pop.Dimen; k++ { + if crossPoints[k] { + trial[k] = currentIndividual.CurX[k] + } + } + + f := bench.Functions[pop.Problem] + currentFitness := f(currentIndividual.CurX) + trialFitness := f(trial) + + // replace if better. + if trialFitness < currentFitness { + pop.Population[i].CurX = trial + } + } +} + +// NewDE returns a pointer to a new, uninitialised DE instance. +func NewDE() *DE { + return &DE{} +} diff --git a/algo/ga/jde.go b/algo/ga/jde.go index 9e846ca..e3b7f45 100644 --- a/algo/ga/jde.go +++ b/algo/ga/jde.go @@ -76,7 +76,14 @@ const ( var jDELogger = newLogger(" *** δ jDE:") // Init initialises the jDE algorithm, performs sanity checks on the inputs. -func (j *JDE) Init(generations, benchMinIters, mutStrategy, adptScheme, np int, f, cr float64, dimensions []int, bench string, ch chan []stats.Stats, chAlgoMeans chan *stats.AlgoBenchMean) { +func (j *JDE) Init( + generations, benchMinIters, mutStrategy, adptScheme, np int, + f, cr float64, + dimensions []int, + bench string, + ch chan []stats.Stats, + chAlgoMeans chan *stats.AlgoBenchMean, +) { if j == nil { jDELogger.Fatalln("jDE needs to be initialised before calling Run(), exiting...") } @@ -124,8 +131,16 @@ func (j *JDE) Init(generations, benchMinIters, mutStrategy, adptScheme, np int, rngsrc := rand.NewSource(uint64(rand.Int63())) - j.rngF = distuv.Uniform{Min: fMin, Max: fMax, Src: rngsrc} - j.rngCR = distuv.Uniform{Min: crMin, Max: crMax, Src: rngsrc} + j.rngF = distuv.Uniform{ + Min: fMin, + Max: fMax, + Src: rngsrc, + } + j.rngCR = distuv.Uniform{ + Min: crMin, + Max: crMax, + Src: rngsrc, + } j.Dimensions = dimensions j.BenchName = bench @@ -139,12 +154,25 @@ func (j *JDE) Init(generations, benchMinIters, mutStrategy, adptScheme, np int, // InitAndRun initialises the jDE algorithm, performs sanity checks on the // inputs and calls the Run method. -func (j *JDE) InitAndRun(generations, benchMinIters, mutStrategy, adptScheme, np int, f, cr float64, dimensions []int, bench string, ch chan []stats.Stats, chAlgoMeans chan *stats.AlgoBenchMean) { +func (j *JDE) InitAndRun( + generations, benchMinIters, mutStrategy, adptScheme, np int, + f, cr float64, + dimensions []int, + bench string, + ch chan []stats.Stats, + chAlgoMeans chan *stats.AlgoBenchMean, +) { if j == nil { jDELogger.Fatalln("jDE is nil, NewjDE() needs to be called first. exiting...") } - j.Init(generations, benchMinIters, mutStrategy, adptScheme, np, f, cr, dimensions, bench, ch, chAlgoMeans) + j.Init( + generations, benchMinIters, mutStrategy, adptScheme, np, + f, cr, + dimensions, + bench, + ch, chAlgoMeans, + ) j.Run() } diff --git a/bench/cec2020/basicFunctions.go b/bench/cec2020/basicFunctions.go index 3880fe1..2deffc0 100644 --- a/bench/cec2020/basicFunctions.go +++ b/bench/cec2020/basicFunctions.go @@ -153,6 +153,7 @@ func SchwefelModified(x []float64) float64 { case zi < -500: // g(zi) + sum += (math.Mod(math.Abs(zi), 500)-500)*math.Sin(math.Sqrt(math.Abs(math.Mod(math.Abs(zi), 500)-500))) - (math.Pow(zi-500, 2) - 10000*fnx) case math.Abs(zi) <= 500: // g(zi) diff --git a/report/comparisonOfMeans.go b/report/comparisonOfMeans.go index fe8ed73..aed0a04 100644 --- a/report/comparisonOfMeans.go +++ b/report/comparisonOfMeans.go @@ -35,9 +35,10 @@ func SaveComparisonOfMeans(p PicList, benchCount int) { // split the slice to smaller, per-bench slices. for i := 0; i < len(p.Pics); i += benchCount { - pL := &PicList{} - pL.Pics = p.Pics[i : i+benchCount] - pL.Bench = p.Pics[i].Bench + pL := &PicList{ + Pics: p.Pics[i : i+benchCount], + Bench: p.Pics[i].Bench, + } benchPicLists = append(benchPicLists, *pL) } diff --git a/report/pics.tmpl b/report/pics.tmpl index 4a650e9..3ee0a2b 100644 --- a/report/pics.tmpl +++ b/report/pics.tmpl @@ -11,7 +11,7 @@ \begin{figure}[h!] \centering {{- range $i, $v := .Pics }} - \begin{subfigure}{0.30\textwidth} + \begin{subfigure}{0.46\textwidth} % note: this accomodates 3 plots a row comfortably..should the requirements % change, this would have to be reworked. % {\includesvg[scale=0.45]{ {{- printf "%s" $v.FilePath -}} }} @@ -24,7 +24,7 @@ {{- end -}} % \newline {{ range $k, $w := .PicsMean }} - \begin{subfigure}{0.30\textwidth} + \begin{subfigure}{0.46\textwidth} \vspace{2em} % {\includesvg[scale=0.45]{ {{- printf "%s" $w.FilePath -}} }} % using .pdf diff --git a/run.go b/run.go index 8747cc5..d9d15d0 100644 --- a/run.go +++ b/run.go @@ -25,7 +25,8 @@ var ( jDE = flag.Bool("jde", false, "run Differential Evolution algorithm with parameter self adaptation") // run CEC2020 jDE by default. - c2jDE = flag.Bool("c2jde", true, "run CEC2020 version of the Differential Evolution algorithm with parameter self adaptation") + c2jDE = flag.Bool("c2jde", true, "run CEC2020 version of the Differential Evolution algorithm with parameter self adaptation") + c2SOMAT3A = flag.Bool("c2somat3a", false, "run CEC2020 version of the SOMA Team-to-Team Adaptive (T3A)") ) func run() { @@ -34,7 +35,7 @@ func run() { flag.Parse() if *generate { - if !*jDE && !*c2jDE && !*sHC && !*rS { + if !*jDE && !*c2jDE && !*c2SOMAT3A && !*sHC && !*rS { log.Println("at least one algo needs to be specified, exiting...") return @@ -56,6 +57,12 @@ func run() { go algo.DoCEC2020jDE(&wg, &m) } + if *c2SOMAT3A { + wg.Add(1) + + go algo.DoCEC2020SOMAT3A(&wg, &m) + } + if *rS { wg.Add(1) @@ -74,9 +81,17 @@ func run() { wg.Wait() - // pL, benchCount := algo.PrepComparisonOfMeans(&wg) + var pL *report.PicList - // report.SaveComparisonOfMeans(*pL, benchCount) + var benchCount int + + if *c2jDE && *c2SOMAT3A { + pL, benchCount = algo.PrepCEC2020ComparisonOfMeans(&wg) + } else { + pL, benchCount = algo.PrepComparisonOfMeans(&wg) + } + + report.SaveComparisonOfMeans(*pL, benchCount) report.SaveTexAllPics() report.SaveTexAllTables() }