From 4daac79f231248ea0fb867bb3538c2979b3292e6 Mon Sep 17 00:00:00 2001 From: Connor Olding Date: Sat, 6 May 2023 20:10:20 -0700 Subject: [PATCH] import `optimizers` from evolopy --- evolopy/BAT.py | 110 ++++++++++++ evolopy/CS.py | 157 +++++++++++++++++ evolopy/DE.py | 124 ++++++++++++++ evolopy/FFA.py | 139 +++++++++++++++ evolopy/GA.py | 439 ++++++++++++++++++++++++++++++++++++++++++++++++ evolopy/GWO.py | 145 ++++++++++++++++ evolopy/HHO.py | 211 +++++++++++++++++++++++ evolopy/JAYA.py | 118 +++++++++++++ evolopy/MFO.py | 179 ++++++++++++++++++++ evolopy/MVO.py | 168 ++++++++++++++++++ evolopy/PSO.py | 103 ++++++++++++ evolopy/SCA.py | 93 ++++++++++ evolopy/SSA.py | 125 ++++++++++++++ evolopy/WOA.py | 132 +++++++++++++++ 14 files changed, 2243 insertions(+) create mode 100644 evolopy/BAT.py create mode 100644 evolopy/CS.py create mode 100644 evolopy/DE.py create mode 100644 evolopy/FFA.py create mode 100644 evolopy/GA.py create mode 100644 evolopy/GWO.py create mode 100644 evolopy/HHO.py create mode 100644 evolopy/JAYA.py create mode 100644 evolopy/MFO.py create mode 100644 evolopy/MVO.py create mode 100644 evolopy/PSO.py create mode 100644 evolopy/SCA.py create mode 100644 evolopy/SSA.py create mode 100644 evolopy/WOA.py diff --git a/evolopy/BAT.py b/evolopy/BAT.py new file mode 100644 index 0000000..e51430f --- /dev/null +++ b/evolopy/BAT.py @@ -0,0 +1,110 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu May 26 02:00:55 2016 + +@author: hossam +""" +import math +import numpy +import random +import time +from solution import solution + + +def BAT(objf, lb, ub, dim, N, Max_iteration): + + n = N + # Population size + + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + N_gen = Max_iteration # Number of generations + + A = 0.5 + # Loudness (constant or decreasing) + r = 0.5 + # Pulse rate (constant or decreasing) + + Qmin = 0 # Frequency minimum + Qmax = 2 # Frequency maximum + + d = dim # Number of dimensions + + # Initializing arrays + Q = numpy.zeros(n) # Frequency + v = numpy.zeros((n, d)) # Velocities + Convergence_curve = [] + + # Initialize the population/solutions + Sol = numpy.zeros((n, d)) + for i in range(dim): + Sol[:, i] = numpy.random.rand(n) * (ub[i] - lb[i]) + lb[i] + + S = numpy.zeros((n, d)) + S = numpy.copy(Sol) + Fitness = numpy.zeros(n) + + # initialize solution for the final results + s = solution() + print('BAT is optimizing "' + objf.__name__ + '"') + + # Initialize timer for the experiment + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + # Evaluate initial random solutions + for i in range(0, n): + Fitness[i] = objf(Sol[i, :]) + + # Find the initial best solution and minimum fitness + I = numpy.argmin(Fitness) + best = Sol[I, :] + fmin = min(Fitness) + + # Main loop + for t in range(0, N_gen): + + # Loop over all bats(solutions) + for i in range(0, n): + Q[i] = Qmin + (Qmin - Qmax) * random.random() + v[i, :] = v[i, :] + (Sol[i, :] - best) * Q[i] + S[i, :] = Sol[i, :] + v[i, :] + + # Check boundaries + for j in range(d): + Sol[i, j] = numpy.clip(Sol[i, j], lb[j], ub[j]) + + # Pulse rate + if random.random() > r: + S[i, :] = best + 0.001 * numpy.random.randn(d) + + # Evaluate new solutions + Fnew = objf(S[i, :]) + + # Update if the solution improves + if (Fnew <= Fitness[i]) and (random.random() < A): + Sol[i, :] = numpy.copy(S[i, :]) + Fitness[i] = Fnew + + # Update the current best solution + if Fnew <= fmin: + best = numpy.copy(S[i, :]) + fmin = Fnew + + # update convergence curve + Convergence_curve.append(fmin) + + if t % 1 == 0: + print(["At iteration " + str(t) + " the best fitness is " + str(fmin)]) + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "BAT" + s.bestIndividual = best + s.objfname = objf.__name__ + + return s diff --git a/evolopy/CS.py b/evolopy/CS.py new file mode 100644 index 0000000..09efe59 --- /dev/null +++ b/evolopy/CS.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue May 24 13:13:28 2016 + +@author: Hossam Faris +""" +import math +import numpy +import random +import time +from solution import solution + + +def get_cuckoos(nest, best, lb, ub, n, dim): + + # perform Levy flights + tempnest = numpy.zeros((n, dim)) + tempnest = numpy.array(nest) + beta = 3 / 2 + sigma = ( + math.gamma(1 + beta) + * math.sin(math.pi * beta / 2) + / (math.gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2)) + ) ** (1 / beta) + + s = numpy.zeros(dim) + for j in range(0, n): + s = nest[j, :] + u = numpy.random.randn(len(s)) * sigma + v = numpy.random.randn(len(s)) + step = u / abs(v) ** (1 / beta) + + stepsize = 0.01 * (step * (s - best)) + + s = s + stepsize * numpy.random.randn(len(s)) + + for k in range(dim): + tempnest[j, k] = numpy.clip(s[k], lb[k], ub[k]) + + return tempnest + + +def get_best_nest(nest, newnest, fitness, n, dim, objf): + # Evaluating all new solutions + tempnest = numpy.zeros((n, dim)) + tempnest = numpy.copy(nest) + + for j in range(0, n): + # for j=1:size(nest,1), + fnew = objf(newnest[j, :]) + if fnew <= fitness[j]: + fitness[j] = fnew + tempnest[j, :] = newnest[j, :] + + # Find the current best + + fmin = min(fitness) + K = numpy.argmin(fitness) + bestlocal = tempnest[K, :] + + return fmin, bestlocal, tempnest, fitness + + +# Replace some nests by constructing new solutions/nests +def empty_nests(nest, pa, n, dim): + + # Discovered or not + tempnest = numpy.zeros((n, dim)) + + K = numpy.random.uniform(0, 1, (n, dim)) > pa + + stepsize = random.random() * ( + nest[numpy.random.permutation(n), :] - nest[numpy.random.permutation(n), :] + ) + + tempnest = nest + stepsize * K + + return tempnest + + +########################################################################## + + +def CS(objf, lb, ub, dim, n, N_IterTotal): + + # lb=-1 + # ub=1 + # n=50 + # N_IterTotal=1000 + # dim=30 + + # Discovery rate of alien eggs/solutions + pa = 0.25 + + nd = dim + + # Lb=[lb]*nd + # Ub=[ub]*nd + convergence = [] + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # RInitialize nests randomely + nest = numpy.zeros((n, dim)) + for i in range(dim): + nest[:, i] = numpy.random.uniform(0, 1, n) * (ub[i] - lb[i]) + lb[i] + + new_nest = numpy.zeros((n, dim)) + new_nest = numpy.copy(nest) + + bestnest = [0] * dim + + fitness = numpy.zeros(n) + fitness.fill(float("inf")) + + s = solution() + + print('CS is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + fmin, bestnest, nest, fitness = get_best_nest(nest, new_nest, fitness, n, dim, objf) + convergence = [] + # Main loop counter + for iter in range(0, N_IterTotal): + # Generate new solutions (but keep the current best) + + new_nest = get_cuckoos(nest, bestnest, lb, ub, n, dim) + + # Evaluate new solutions and find best + fnew, best, nest, fitness = get_best_nest(nest, new_nest, fitness, n, dim, objf) + + new_nest = empty_nests(new_nest, pa, n, dim) + + # Evaluate new solutions and find best + fnew, best, nest, fitness = get_best_nest(nest, new_nest, fitness, n, dim, objf) + + if fnew < fmin: + fmin = fnew + bestnest = best + + if iter % 10 == 0: + print(["At iteration " + str(iter) + " the best fitness is " + str(fmin)]) + convergence.append(fmin) + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence + s.optimizer = "CS" + s.bestIndividual = bestnest + s.objfname = objf.__name__ + + return s diff --git a/evolopy/DE.py b/evolopy/DE.py new file mode 100644 index 0000000..bc206e2 --- /dev/null +++ b/evolopy/DE.py @@ -0,0 +1,124 @@ +import random +import numpy +import time +from solution import solution + + +# Differential Evolution (DE) +# mutation factor = [0.5, 2] +# crossover_ratio = [0,1] +def DE(objf, lb, ub, dim, PopSize, iters): + + mutation_factor = 0.5 + crossover_ratio = 0.7 + stopping_func = None + + # convert lb, ub to array + if not isinstance(lb, list): + lb = [lb for _ in range(dim)] + ub = [ub for _ in range(dim)] + + # solution + s = solution() + + s.best = float("inf") + + # initialize population + population = [] + + population_fitness = numpy.array([float("inf") for _ in range(PopSize)]) + + for p in range(PopSize): + sol = [] + for d in range(dim): + d_val = random.uniform(lb[d], ub[d]) + sol.append(d_val) + + population.append(sol) + + population = numpy.array(population) + + # calculate fitness for all the population + for i in range(PopSize): + fitness = objf(population[i, :]) + population_fitness[p] = fitness + # s.func_evals += 1 + + # is leader ? + if fitness < s.best: + s.best = fitness + s.leader_solution = population[i, :] + + convergence_curve = numpy.zeros(iters) + # start work + print('DE is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + t = 0 + while t < iters: + # should i stop + if stopping_func is not None and stopping_func(s.best, s.leader_solution, t): + break + + # loop through population + for i in range(PopSize): + # 1. Mutation + + # select 3 random solution except current solution + ids_except_current = [_ for _ in range(PopSize) if _ != i] + id_1, id_2, id_3 = random.sample(ids_except_current, 3) + + mutant_sol = [] + for d in range(dim): + d_val = population[id_1, d] + mutation_factor * ( + population[id_2, d] - population[id_3, d] + ) + + # 2. Recombination + rn = random.uniform(0, 1) + if rn > crossover_ratio: + d_val = population[i, d] + + # add dimension value to the mutant solution + mutant_sol.append(d_val) + + # 3. Replacement / Evaluation + + # clip new solution (mutant) + mutant_sol = numpy.clip(mutant_sol, lb, ub) + + # calc fitness + mutant_fitness = objf(mutant_sol) + # s.func_evals += 1 + + # replace if mutant_fitness is better + if mutant_fitness < population_fitness[i]: + population[i, :] = mutant_sol + population_fitness[i] = mutant_fitness + + # update leader + if mutant_fitness < s.best: + s.best = mutant_fitness + s.leader_solution = mutant_sol + + convergence_curve[t] = s.best + if t % 1 == 0: + print( + ["At iteration " + str(t + 1) + " the best fitness is " + str(s.best)] + ) + + # increase iterations + t = t + 1 + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence_curve + s.optimizer = "DE" + s.bestIndividual = s.leader_solution + s.objfname = objf.__name__ + + # return solution + return s diff --git a/evolopy/FFA.py b/evolopy/FFA.py new file mode 100644 index 0000000..3a87ac8 --- /dev/null +++ b/evolopy/FFA.py @@ -0,0 +1,139 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun May 29 00:49:35 2016 + +@author: hossam +""" + +#% ======================================================== % +#% Files of the Matlab programs included in the book: % +#% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms, % +#% Second Edition, Luniver Press, (2010). www.luniver.com % +#% ======================================================== % +# +#% -------------------------------------------------------- % +#% Firefly Algorithm for constrained optimization using % +#% for the design of a spring (benchmark) % +#% by Xin-She Yang (Cambridge University) Copyright @2009 % +#% -------------------------------------------------------- % + +import numpy +import math +import time +from solution import solution + + +def alpha_new(alpha, NGen): + #% alpha_n=alpha_0(1-delta)^NGen=10^(-4); + #% alpha_0=0.9 + delta = 1 - (10 ** (-4) / 0.9) ** (1 / NGen) + alpha = (1 - delta) * alpha + return alpha + + +def FFA(objf, lb, ub, dim, n, MaxGeneration): + + # General parameters + + # n=50 #number of fireflies + # dim=30 #dim + # lb=-50 + # ub=50 + # MaxGeneration=500 + + # FFA parameters + alpha = 0.5 # Randomness 0--1 (highly random) + betamin = 0.20 # minimum value of beta + gamma = 1 # Absorption coefficient + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + zn = numpy.ones(n) + zn.fill(float("inf")) + + # ns(i,:)=Lb+(Ub-Lb).*rand(1,d); + ns = numpy.zeros((n, dim)) + for i in range(dim): + ns[:, i] = numpy.random.uniform(0, 1, n) * (ub[i] - lb[i]) + lb[i] + Lightn = numpy.ones(n) + Lightn.fill(float("inf")) + + # [ns,Lightn]=init_ffa(n,d,Lb,Ub,u0) + + convergence = [] + s = solution() + + print('FFA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + # Main loop + for k in range(0, MaxGeneration): # start iterations + + #% This line of reducing alpha is optional + alpha = alpha_new(alpha, MaxGeneration) + + #% Evaluate new solutions (for all n fireflies) + for i in range(0, n): + zn[i] = objf(ns[i, :]) + Lightn[i] = zn[i] + + # Ranking fireflies by their light intensity/objectives + + Lightn = numpy.sort(zn) + Index = numpy.argsort(zn) + ns = ns[Index, :] + + # Find the current best + nso = ns + Lighto = Lightn + nbest = ns[0, :] + Lightbest = Lightn[0] + + #% For output only + fbest = Lightbest + + #% Move all fireflies to the better locations + # [ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,nbest,... + # Lightbest,alpha,betamin,gamma,Lb,Ub); + scale = [] + for b in range(dim): + scale.append(abs(ub[b] - lb[b])) + scale = numpy.array(scale) + for i in range(0, n): + # The attractiveness parameter beta=exp(-gamma*r) + for j in range(0, n): + r = numpy.sqrt(numpy.sum((ns[i, :] - ns[j, :]) ** 2)) + # r=1 + # Update moves + if Lightn[i] > Lighto[j]: # Brighter and more attractive + beta0 = 1 + beta = (beta0 - betamin) * math.exp(-gamma * r ** 2) + betamin + tmpf = alpha * (numpy.random.rand(dim) - 0.5) * scale + ns[i, :] = ns[i, :] * (1 - beta) + nso[j, :] * beta + tmpf + + # ns=numpy.clip(ns, lb, ub) + + convergence.append(fbest) + + IterationNumber = k + BestQuality = fbest + + if k % 1 == 0: + print( + ["At iteration " + str(k) + " the best fitness is " + str(BestQuality)] + ) + # + ####################### End main loop + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence + s.optimizer = "FFA" + s.bestIndividual = nbest + s.objfname = objf.__name__ + + return s diff --git a/evolopy/GA.py b/evolopy/GA.py new file mode 100644 index 0000000..485b140 --- /dev/null +++ b/evolopy/GA.py @@ -0,0 +1,439 @@ +""" +Created on Sat Feb 24 20:18:05 2019 + +@author: Raneem +""" +import numpy +import random +import time +import sys + +from solution import solution + + +def crossoverPopulaton(population, scores, popSize, crossoverProbability, keep): + """ + The crossover of all individuals + + Parameters + ---------- + population : list + The list of individuals + scores : list + The list of fitness values for each individual + popSize: int + Number of chrmosome in a population + crossoverProbability: float + The probability of crossing a pair of individuals + keep: int + Number of best individuals to keep without mutating for the next generation + + + Returns + ------- + N/A + """ + # initialize a new population + newPopulation = numpy.empty_like(population) + newPopulation[0:keep] = population[0:keep] + # Create pairs of parents. The number of pairs equals the number of individuals divided by 2 + for i in range(keep, popSize, 2): + # pair of parents selection + parent1, parent2 = pairSelection(population, scores, popSize) + crossoverLength = min(len(parent1), len(parent2)) + parentsCrossoverProbability = random.uniform(0.0, 1.0) + if parentsCrossoverProbability < crossoverProbability: + offspring1, offspring2 = crossover(crossoverLength, parent1, parent2) + else: + offspring1 = parent1.copy() + offspring2 = parent2.copy() + + # Add offsprings to population + newPopulation[i] = numpy.copy(offspring1) + newPopulation[i + 1] = numpy.copy(offspring2) + + return newPopulation + + +def mutatePopulaton(population, popSize, mutationProbability, keep, lb, ub): + """ + The mutation of all individuals + + Parameters + ---------- + population : list + The list of individuals + popSize: int + Number of chrmosome in a population + mutationProbability: float + The probability of mutating an individual + keep: int + Number of best individuals to keep without mutating for the next generation + lb: list + lower bound limit list + ub: list + Upper bound limit list + + Returns + ------- + N/A + """ + for i in range(keep, popSize): + # Mutation + offspringMutationProbability = random.uniform(0.0, 1.0) + if offspringMutationProbability < mutationProbability: + mutation(population[i], len(population[i]), lb, ub) + + +def elitism(population, scores, bestIndividual, bestScore): + """ + This melitism operator of the population + + Parameters + ---------- + population : list + The list of individuals + scores : list + The list of fitness values for each individual + bestIndividual : list + An individual of the previous generation having the best fitness value + bestScore : float + The best fitness value of the previous generation + + Returns + ------- + N/A + """ + + # get the worst individual + worstFitnessId = selectWorstIndividual(scores) + + # replace worst cromosome with best one from previous generation if its fitness is less than the other + if scores[worstFitnessId] > bestScore: + population[worstFitnessId] = numpy.copy(bestIndividual) + scores[worstFitnessId] = numpy.copy(bestScore) + + +def selectWorstIndividual(scores): + """ + It is used to get the worst individual in a population based n the fitness value + + Parameters + ---------- + scores : list + The list of fitness values for each individual + + Returns + ------- + int + maxFitnessId: The individual id of the worst fitness value + """ + + maxFitnessId = numpy.where(scores == numpy.max(scores)) + maxFitnessId = maxFitnessId[0][0] + return maxFitnessId + + +def pairSelection(population, scores, popSize): + """ + This is used to select one pair of parents using roulette Wheel Selection mechanism + + Parameters + ---------- + population : list + The list of individuals + scores : list + The list of fitness values for each individual + popSize: int + Number of chrmosome in a population + + Returns + ------- + list + parent1: The first parent individual of the pair + list + parent2: The second parent individual of the pair + """ + parent1Id = rouletteWheelSelectionId(scores, popSize) + parent1 = population[parent1Id].copy() + + parent2Id = rouletteWheelSelectionId(scores, popSize) + parent2 = population[parent2Id].copy() + + return parent1, parent2 + + +def rouletteWheelSelectionId(scores, popSize): + """ + A roulette Wheel Selection mechanism for selecting an individual + + Parameters + ---------- + scores : list + The list of fitness values for each individual + popSize: int + Number of chrmosome in a population + + Returns + ------- + id + individualId: The id of the individual selected + """ + + ##reverse score because minimum value should have more chance of selection + reverse = max(scores) + min(scores) + reverseScores = reverse - scores.copy() + sumScores = sum(reverseScores) + pick = random.uniform(0, sumScores) + current = 0 + for individualId in range(popSize): + current += reverseScores[individualId] + if current > pick: + return individualId + + +def crossover(individualLength, parent1, parent2): + """ + The crossover operator of a two individuals + + Parameters + ---------- + individualLength: int + The maximum index of the crossover + parent1 : list + The first parent individual of the pair + parent2 : list + The second parent individual of the pair + + Returns + ------- + list + offspring1: The first updated parent individual of the pair + list + offspring2: The second updated parent individual of the pair + """ + + # The point at which crossover takes place between two parents. + crossover_point = random.randint(0, individualLength - 1) + # The new offspring will have its first half of its genes taken from the first parent and second half of its genes taken from the second parent. + offspring1 = numpy.concatenate( + [parent1[0:crossover_point], parent2[crossover_point:]] + ) + # The new offspring will have its first half of its genes taken from the second parent and second half of its genes taken from the first parent. + offspring2 = numpy.concatenate( + [parent2[0:crossover_point], parent1[crossover_point:]] + ) + + return offspring1, offspring2 + + +def mutation(offspring, individualLength, lb, ub): + """ + The mutation operator of a single individual + + Parameters + ---------- + offspring : list + A generated individual after the crossover + individualLength: int + The maximum index of the crossover + lb: list + lower bound limit list + ub: list + Upper bound limit list + + Returns + ------- + N/A + """ + mutationIndex = random.randint(0, individualLength - 1) + mutationValue = random.uniform(lb[mutationIndex], ub[mutationIndex]) + offspring[mutationIndex] = mutationValue + + +def clearDups(Population, lb, ub): + + """ + It removes individuals duplicates and replace them with random ones + + Parameters + ---------- + objf : function + The objective function selected + lb: list + lower bound limit list + ub: list + Upper bound limit list + + Returns + ------- + list + newPopulation: the updated list of individuals + """ + newPopulation = numpy.unique(Population, axis=0) + oldLen = len(Population) + newLen = len(newPopulation) + if newLen < oldLen: + nDuplicates = oldLen - newLen + newPopulation = numpy.append( + newPopulation, + numpy.random.uniform(0, 1, (nDuplicates, len(Population[0]))) + * (numpy.array(ub) - numpy.array(lb)) + + numpy.array(lb), + axis=0, + ) + + return newPopulation + + +def calculateCost(objf, population, popSize, lb, ub): + + """ + It calculates the fitness value of each individual in the population + + Parameters + ---------- + objf : function + The objective function selected + population : list + The list of individuals + popSize: int + Number of chrmosomes in a population + lb: list + lower bound limit list + ub: list + Upper bound limit list + + Returns + ------- + list + scores: fitness values of all individuals in the population + """ + scores = numpy.full(popSize, numpy.inf) + + # Loop through individuals in population + for i in range(0, popSize): + # Return back the search agents that go beyond the boundaries of the search space + population[i] = numpy.clip(population[i], lb, ub) + + # Calculate objective function for each search agent + scores[i] = objf(population[i, :]) + + return scores + + +def sortPopulation(population, scores): + """ + This is used to sort the population according to the fitness values of the individuals + + Parameters + ---------- + population : list + The list of individuals + scores : list + The list of fitness values for each individual + + Returns + ------- + list + population: The new sorted list of individuals + list + scores: The new sorted list of fitness values of the individuals + """ + sortedIndices = scores.argsort() + population = population[sortedIndices] + scores = scores[sortedIndices] + + return population, scores + + +def GA(objf, lb, ub, dim, popSize, iters): + + """ + This is the main method which implements GA + + Parameters + ---------- + objf : function + The objective function selected + lb: list + lower bound limit list + ub: list + Upper bound limit list + dim: int + The dimension of the indivisual + popSize: int + Number of chrmosomes in a population + iters: int + Number of iterations / generations of GA + + Returns + ------- + obj + s: The solution obtained from running the algorithm + """ + + cp = 1 # crossover Probability + mp = 0.01 # Mutation Probability + keep = 2 + # elitism parameter: how many of the best individuals to keep from one generation to the next + + s = solution() + + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + bestIndividual = numpy.zeros(dim) + scores = numpy.random.uniform(0.0, 1.0, popSize) + bestScore = float("inf") + + ga = numpy.zeros((popSize, dim)) + for i in range(dim): + ga[:, i] = numpy.random.uniform(0, 1, popSize) * (ub[i] - lb[i]) + lb[i] + convergence_curve = numpy.zeros(iters) + + print('GA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + for l in range(iters): + + # crossover + ga = crossoverPopulaton(ga, scores, popSize, cp, keep) + + # mutation + mutatePopulaton(ga, popSize, mp, keep, lb, ub) + + ga = clearDups(ga, lb, ub) + + scores = calculateCost(objf, ga, popSize, lb, ub) + + bestScore = min(scores) + + # Sort from best to worst + ga, scores = sortPopulation(ga, scores) + + convergence_curve[l] = bestScore + + if l % 1 == 0: + print( + [ + "At iteration " + + str(l + 1) + + " the best fitness is " + + str(bestScore) + ] + ) + + timerEnd = time.time() + s.bestIndividual = bestIndividual + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence_curve + s.optimizer = "GA" + s.objfname = objf.__name__ + + return s diff --git a/evolopy/GWO.py b/evolopy/GWO.py new file mode 100644 index 0000000..4aee47b --- /dev/null +++ b/evolopy/GWO.py @@ -0,0 +1,145 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon May 16 00:27:50 2016 + +@author: Hossam Faris +""" + +import random +import numpy +import math +from solution import solution +import time + + +def GWO(objf, lb, ub, dim, SearchAgents_no, Max_iter): + + # Max_iter=1000 + # lb=-100 + # ub=100 + # dim=30 + # SearchAgents_no=5 + + # initialize alpha, beta, and delta_pos + Alpha_pos = numpy.zeros(dim) + Alpha_score = float("inf") + + Beta_pos = numpy.zeros(dim) + Beta_score = float("inf") + + Delta_pos = numpy.zeros(dim) + Delta_score = float("inf") + + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # Initialize the positions of search agents + Positions = numpy.zeros((SearchAgents_no, dim)) + for i in range(dim): + Positions[:, i] = ( + numpy.random.uniform(0, 1, SearchAgents_no) * (ub[i] - lb[i]) + lb[i] + ) + + Convergence_curve = numpy.zeros(Max_iter) + s = solution() + + # Loop counter + print('GWO is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + # Main loop + for l in range(0, Max_iter): + for i in range(0, SearchAgents_no): + + # Return back the search agents that go beyond the boundaries of the search space + for j in range(dim): + Positions[i, j] = numpy.clip(Positions[i, j], lb[j], ub[j]) + + # Calculate objective function for each search agent + fitness = objf(Positions[i, :]) + + # Update Alpha, Beta, and Delta + if fitness < Alpha_score: + Delta_score = Beta_score # Update delte + Delta_pos = Beta_pos.copy() + Beta_score = Alpha_score # Update beta + Beta_pos = Alpha_pos.copy() + Alpha_score = fitness + # Update alpha + Alpha_pos = Positions[i, :].copy() + + if fitness > Alpha_score and fitness < Beta_score: + Delta_score = Beta_score # Update delte + Delta_pos = Beta_pos.copy() + Beta_score = fitness # Update beta + Beta_pos = Positions[i, :].copy() + + if fitness > Alpha_score and fitness > Beta_score and fitness < Delta_score: + Delta_score = fitness # Update delta + Delta_pos = Positions[i, :].copy() + + a = 2 - l * ((2) / Max_iter) + # a decreases linearly fron 2 to 0 + + # Update the Position of search agents including omegas + for i in range(0, SearchAgents_no): + for j in range(0, dim): + + r1 = random.random() # r1 is a random number in [0,1] + r2 = random.random() # r2 is a random number in [0,1] + + A1 = 2 * a * r1 - a + # Equation (3.3) + C1 = 2 * r2 + # Equation (3.4) + + D_alpha = abs(C1 * Alpha_pos[j] - Positions[i, j]) + # Equation (3.5)-part 1 + X1 = Alpha_pos[j] - A1 * D_alpha + # Equation (3.6)-part 1 + + r1 = random.random() + r2 = random.random() + + A2 = 2 * a * r1 - a + # Equation (3.3) + C2 = 2 * r2 + # Equation (3.4) + + D_beta = abs(C2 * Beta_pos[j] - Positions[i, j]) + # Equation (3.5)-part 2 + X2 = Beta_pos[j] - A2 * D_beta + # Equation (3.6)-part 2 + + r1 = random.random() + r2 = random.random() + + A3 = 2 * a * r1 - a + # Equation (3.3) + C3 = 2 * r2 + # Equation (3.4) + + D_delta = abs(C3 * Delta_pos[j] - Positions[i, j]) + # Equation (3.5)-part 3 + X3 = Delta_pos[j] - A3 * D_delta + # Equation (3.5)-part 3 + + Positions[i, j] = (X1 + X2 + X3) / 3 # Equation (3.7) + + Convergence_curve[l] = Alpha_score + + if l % 1 == 0: + print(["At iteration " + str(l) + " the best fitness is " + str(Alpha_score)]) + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "GWO" + s.bestIndividual = Alpha_pos + s.objfname = objf.__name__ + + return s diff --git a/evolopy/HHO.py b/evolopy/HHO.py new file mode 100644 index 0000000..9b37af8 --- /dev/null +++ b/evolopy/HHO.py @@ -0,0 +1,211 @@ +# -*- coding: utf-8 -*- +""" +Created on Thirsday March 21 2019 + +@author: +% _____________________________________________________ +% Main paper: +% Harris hawks optimization: Algorithm and applications +% Ali Asghar Heidari, Seyedali Mirjalili, Hossam Faris, Ibrahim Aljarah, Majdi Mafarja, Huiling Chen +% Future Generation Computer Systems, +% DOI: https://doi.org/10.1016/j.future.2019.02.028 +% _____________________________________________________ + +""" +import random +import numpy +import math +from solution import solution +import time + + +def HHO(objf, lb, ub, dim, SearchAgents_no, Max_iter): + + # dim=30 + # SearchAgents_no=50 + # lb=-100 + # ub=100 + # Max_iter=500 + + # initialize the location and Energy of the rabbit + Rabbit_Location = numpy.zeros(dim) + Rabbit_Energy = float("inf") # change this to -inf for maximization problems + + if not isinstance(lb, list): + lb = [lb for _ in range(dim)] + ub = [ub for _ in range(dim)] + lb = numpy.asarray(lb) + ub = numpy.asarray(ub) + + # Initialize the locations of Harris' hawks + X = numpy.asarray( + [x * (ub - lb) + lb for x in numpy.random.uniform(0, 1, (SearchAgents_no, dim))] + ) + + # Initialize convergence + convergence_curve = numpy.zeros(Max_iter) + + ############################ + s = solution() + + print('HHO is now tackling "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + ############################ + + t = 0 # Loop counter + + # Main loop + while t < Max_iter: + for i in range(0, SearchAgents_no): + + # Check boundries + + X[i, :] = numpy.clip(X[i, :], lb, ub) + + # fitness of locations + fitness = objf(X[i, :]) + + # Update the location of Rabbit + if fitness < Rabbit_Energy: # Change this to > for maximization problem + Rabbit_Energy = fitness + Rabbit_Location = X[i, :].copy() + + E1 = 2 * (1 - (t / Max_iter)) # factor to show the decreaing energy of rabbit + + # Update the location of Harris' hawks + for i in range(0, SearchAgents_no): + + E0 = 2 * random.random() - 1 # -1= 1: + # Harris' hawks perch randomly based on 2 strategy: + q = random.random() + rand_Hawk_index = math.floor(SearchAgents_no * random.random()) + X_rand = X[rand_Hawk_index, :] + if q < 0.5: + # perch based on other family members + X[i, :] = X_rand - random.random() * abs( + X_rand - 2 * random.random() * X[i, :] + ) + + elif q >= 0.5: + # perch on a random tall tree (random site inside group's home range) + X[i, :] = (Rabbit_Location - X.mean(0)) - random.random() * ( + (ub - lb) * random.random() + lb + ) + + # -------- Exploitation phase ------------------- + elif abs(Escaping_Energy) < 1: + # Attacking the rabbit using 4 strategies regarding the behavior of the rabbit + + # phase 1: ----- surprise pounce (seven kills) ---------- + # surprise pounce (seven kills): multiple, short rapid dives by different hawks + + r = random.random() # probablity of each event + + if ( + r >= 0.5 and abs(Escaping_Energy) < 0.5 + ): # Hard besiege Eq. (6) in paper + X[i, :] = (Rabbit_Location) - Escaping_Energy * abs( + Rabbit_Location - X[i, :] + ) + + if ( + r >= 0.5 and abs(Escaping_Energy) >= 0.5 + ): # Soft besiege Eq. (4) in paper + Jump_strength = 2 * ( + 1 - random.random() + ) # random jump strength of the rabbit + X[i, :] = (Rabbit_Location - X[i, :]) - Escaping_Energy * abs( + Jump_strength * Rabbit_Location - X[i, :] + ) + + # phase 2: --------performing team rapid dives (leapfrog movements)---------- + + if ( + r < 0.5 and abs(Escaping_Energy) >= 0.5 + ): # Soft besiege Eq. (10) in paper + # rabbit try to escape by many zigzag deceptive motions + Jump_strength = 2 * (1 - random.random()) + X1 = Rabbit_Location - Escaping_Energy * abs( + Jump_strength * Rabbit_Location - X[i, :] + ) + X1 = numpy.clip(X1, lb, ub) + + if objf(X1) < fitness: # improved move? + X[i, :] = X1.copy() + else: # hawks perform levy-based short rapid dives around the rabbit + X2 = ( + Rabbit_Location + - Escaping_Energy + * abs(Jump_strength * Rabbit_Location - X[i, :]) + + numpy.multiply(numpy.random.randn(dim), Levy(dim)) + ) + X2 = numpy.clip(X2, lb, ub) + if objf(X2) < fitness: + X[i, :] = X2.copy() + if ( + r < 0.5 and abs(Escaping_Energy) < 0.5 + ): # Hard besiege Eq. (11) in paper + Jump_strength = 2 * (1 - random.random()) + X1 = Rabbit_Location - Escaping_Energy * abs( + Jump_strength * Rabbit_Location - X.mean(0) + ) + X1 = numpy.clip(X1, lb, ub) + + if objf(X1) < fitness: # improved move? + X[i, :] = X1.copy() + else: # Perform levy-based short rapid dives around the rabbit + X2 = ( + Rabbit_Location + - Escaping_Energy + * abs(Jump_strength * Rabbit_Location - X.mean(0)) + + numpy.multiply(numpy.random.randn(dim), Levy(dim)) + ) + X2 = numpy.clip(X2, lb, ub) + if objf(X2) < fitness: + X[i, :] = X2.copy() + + convergence_curve[t] = Rabbit_Energy + if t % 1 == 0: + print( + [ + "At iteration " + + str(t) + + " the best fitness is " + + str(Rabbit_Energy) + ] + ) + t = t + 1 + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence_curve + s.optimizer = "HHO" + s.objfname = objf.__name__ + s.best = Rabbit_Energy + s.bestIndividual = Rabbit_Location + + return s + + +def Levy(dim): + beta = 1.5 + sigma = ( + math.gamma(1 + beta) + * math.sin(math.pi * beta / 2) + / (math.gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2)) + ) ** (1 / beta) + u = 0.01 * numpy.random.randn(dim) * sigma + v = numpy.random.randn(dim) + zz = numpy.power(numpy.absolute(v), (1 / beta)) + step = numpy.divide(u, zz) + return step diff --git a/evolopy/JAYA.py b/evolopy/JAYA.py new file mode 100644 index 0000000..84c0536 --- /dev/null +++ b/evolopy/JAYA.py @@ -0,0 +1,118 @@ +""" JAYA Algorithm """ + +import random +import numpy +import math +from solution import solution +import time + + +def JAYA(objf, lb, ub, dim, SearchAgents_no, Max_iter): + + # Best and Worst position initialization + Best_pos = numpy.zeros(dim) + Best_score = float("inf") + + Worst_pos = numpy.zeros(dim) + Worst_score = float(0) + + fitness_matrix = numpy.zeros((SearchAgents_no)) + + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # Initialize the positions of search agents + Positions = numpy.zeros((SearchAgents_no, dim)) + for i in range(dim): + Positions[:, i] = ( + numpy.random.uniform(0, 1, SearchAgents_no) * (ub[i] - lb[i]) + lb[i] + ) + + for i in range(0, SearchAgents_no): + + # Return back the search agents that go beyond the boundaries of the search space + for j in range(dim): + Positions[i, j] = numpy.clip(Positions[i, j], lb[j], ub[j]) + + # Calculate objective function for each search agent + fitness = objf(Positions[i]) + fitness_matrix[i] = fitness + + if fitness < Best_score: + Best_score = fitness # Update Best_Score + Best_pos = Positions[i] + + if fitness > Worst_score: + Worst_score = fitness # Update Worst_Score + Worst_pos = Positions[i] + + Convergence_curve = numpy.zeros(Max_iter) + s = solution() + + # Loop counter + print('JAYA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + # Main loop + for l in range(0, Max_iter): + + # Update the Position of search agents + for i in range(0, SearchAgents_no): + New_Position = numpy.zeros(dim) + for j in range(0, dim): + + # Update r1, r2 + r1 = random.random() + r2 = random.random() + + # JAYA Equation + New_Position[j] = ( + Positions[i][j] + + r1 * (Best_pos[j] - abs(Positions[i, j])) + - r2 * (Worst_pos[j] - abs(Positions[i, j])) + ) + + # checking if New_Position[j] lies in search space + if New_Position[j] > ub[j]: + New_Position[j] = ub[j] + if New_Position[j] < lb[j]: + New_Position[j] = lb[j] + + new_fitness = objf(New_Position) + current_fit = fitness_matrix[i] + + # replacing current element with new element if it has better fitness + if new_fitness < current_fit: + Positions[i] = New_Position + fitness_matrix[i] = new_fitness + + # finding the best and worst element + for i in range(SearchAgents_no): + if fitness_matrix[i] < Best_score: + Best_score = fitness_matrix[i] + Best_pos = Positions[i, :].copy() + + if fitness_matrix[i] > Worst_score: + Worst_score = fitness_matrix[i] + Worst_pos = Positions[i, :].copy() + + Convergence_curve[l] = Best_score + + if l % 1 == 0: + print( + ["At iteration " + str(l) + " the best fitness is " + str(Best_score)] + ) + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "JAYA" + s.bestIndividual = Best_pos + s.objfname = objf.__name__ + + return s diff --git a/evolopy/MFO.py b/evolopy/MFO.py new file mode 100644 index 0000000..4f8248e --- /dev/null +++ b/evolopy/MFO.py @@ -0,0 +1,179 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon May 16 10:42:18 2016 + +@author: hossam +""" + +import random +import numpy +import math +from solution import solution +import time + + +def MFO(objf, lb, ub, dim, N, Max_iteration): + + # Max_iteration=1000 + # lb=-100 + # ub=100 + # dim=30 + #N = 50 # Number of search agents + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # Initialize the positions of moths + Moth_pos = numpy.zeros((N, dim)) + for i in range(dim): + Moth_pos[:, i] = numpy.random.uniform(0, 1, N) * (ub[i] - lb[i]) + lb[i] + Moth_fitness = numpy.full(N, float("inf")) + # Moth_fitness=numpy.fell(float("inf")) + + Convergence_curve = numpy.zeros(Max_iteration) + + sorted_population = numpy.copy(Moth_pos) + fitness_sorted = numpy.zeros(N) + ##################### + best_flames = numpy.copy(Moth_pos) + best_flame_fitness = numpy.zeros(N) + #################### + double_population = numpy.zeros((2 * N, dim)) + double_fitness = numpy.zeros(2 * N) + + double_sorted_population = numpy.zeros((2 * N, dim)) + double_fitness_sorted = numpy.zeros(2 * N) + ######################### + previous_population = numpy.zeros((N, dim)) + previous_fitness = numpy.zeros(N) + + s = solution() + + print('MFO is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + Iteration = 1 + + # Main loop + while Iteration < Max_iteration: + + # Number of flames Eq. (3.14) in the paper + Flame_no = round(N - Iteration * ((N - 1) / Max_iteration)) + + for i in range(0, N): + + # Check if moths go out of the search spaceand bring it back + for j in range(dim): + Moth_pos[i, j] = numpy.clip(Moth_pos[i, j], lb[j], ub[j]) + + # evaluate moths + Moth_fitness[i] = objf(Moth_pos[i, :]) + + if Iteration == 1: + # Sort the first population of moths + fitness_sorted = numpy.sort(Moth_fitness) + I = numpy.argsort(Moth_fitness) + + sorted_population = Moth_pos[I, :] + + # Update the flames + best_flames = sorted_population + best_flame_fitness = fitness_sorted + else: + # + # # Sort the moths + double_population = numpy.concatenate( + (previous_population, best_flames), axis=0 + ) + double_fitness = numpy.concatenate( + (previous_fitness, best_flame_fitness), axis=0 + ) + # + double_fitness_sorted = numpy.sort(double_fitness) + I2 = numpy.argsort(double_fitness) + # + # + for newindex in range(0, 2 * N): + double_sorted_population[newindex, :] = numpy.array( + double_population[I2[newindex], :] + ) + + fitness_sorted = double_fitness_sorted[0:N] + sorted_population = double_sorted_population[0:N, :] + # + # # Update the flames + best_flames = sorted_population + best_flame_fitness = fitness_sorted + + # + # # Update the position best flame obtained so far + Best_flame_score = fitness_sorted[0] + Best_flame_pos = sorted_population[0, :] + # + previous_population = Moth_pos + previous_fitness = Moth_fitness + # + # a linearly dicreases from -1 to -2 to calculate t in Eq. (3.12) + a = -1 + Iteration * ((-1) / Max_iteration) + + # Loop counter + for i in range(0, N): + # + for j in range(0, dim): + if ( + i <= Flame_no + ): # Update the position of the moth with respect to its corresponsing flame + # + # D in Eq. (3.13) + distance_to_flame = abs(sorted_population[i, j] - Moth_pos[i, j]) + b = 1 + t = (a - 1) * random.random() + 1 + # + # % Eq. (3.12) + Moth_pos[i, j] = ( + distance_to_flame * math.exp(b * t) * math.cos(t * 2 * math.pi) + + sorted_population[i, j] + ) + # end + # + if ( + i > Flame_no + ): # Upaate the position of the moth with respct to one flame + # + # % Eq. (3.13) + distance_to_flame = abs(sorted_population[i, j] - Moth_pos[i, j]) + b = 1 + t = (a - 1) * random.random() + 1 + # + # % Eq. (3.12) + Moth_pos[i, j] = ( + distance_to_flame * math.exp(b * t) * math.cos(t * 2 * math.pi) + + sorted_population[Flame_no, j] + ) + + Convergence_curve[Iteration] = Best_flame_score + # Display best fitness along the iteration + if Iteration % 1 == 0: + print( + [ + "At iteration " + + str(Iteration) + + " the best fitness is " + + str(Best_flame_score) + ] + ) + + Iteration = Iteration + 1 + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "MFO" + s.bestIndividual = Best_flame_pos + s.objfname = objf.__name__ + + return s diff --git a/evolopy/MVO.py b/evolopy/MVO.py new file mode 100644 index 0000000..ff2be99 --- /dev/null +++ b/evolopy/MVO.py @@ -0,0 +1,168 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed May 11 17:06:34 2016 + +@author: hossam +""" +import random +import numpy +import time +import math +import sklearn +from numpy import asarray +from sklearn.preprocessing import normalize +from solution import solution + + +def normr(Mat): + """normalize the columns of the matrix + B= normr(A) normalizes the row + the dtype of A is float""" + Mat = Mat.reshape(1, -1) + # Enforce dtype float + if Mat.dtype != "float": + Mat = asarray(Mat, dtype=float) + + # if statement to enforce dtype float + B = normalize(Mat, norm="l2", axis=1) + B = numpy.reshape(B, -1) + return B + + +def randk(t): + if (t % 2) == 0: + s = 0.25 + else: + s = 0.75 + return s + + +def RouletteWheelSelection(weights): + accumulation = numpy.cumsum(weights) + p = random.random() * accumulation[-1] + chosen_index = -1 + for index in range(0, len(accumulation)): + if accumulation[index] > p: + chosen_index = index + break + + choice = chosen_index + + return choice + + +def MVO(objf, lb, ub, dim, N, Max_time): + + "parameters" + # dim=30 + # lb=-100 + # ub=100 + WEP_Max = 1 + WEP_Min = 0.2 + # Max_time=1000 + # N=50 + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + Universes = numpy.zeros((N, dim)) + for i in range(dim): + Universes[:, i] = numpy.random.uniform(0, 1, N) * (ub[i] - lb[i]) + lb[i] + + Sorted_universes = numpy.copy(Universes) + + convergence = numpy.zeros(Max_time) + + Best_universe = [0] * dim + Best_universe_Inflation_rate = float("inf") + + s = solution() + + Time = 1 + ############################################ + print('MVO is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + while Time < Max_time + 1: + + "Eq. (3.3) in the paper" + WEP = WEP_Min + Time * ((WEP_Max - WEP_Min) / Max_time) + + TDR = 1 - (math.pow(Time, 1 / 6) / math.pow(Max_time, 1 / 6)) + + Inflation_rates = [0] * len(Universes) + + for i in range(0, N): + for j in range(dim): + Universes[i, j] = numpy.clip(Universes[i, j], lb[j], ub[j]) + + Inflation_rates[i] = objf(Universes[i, :]) + + if Inflation_rates[i] < Best_universe_Inflation_rate: + + Best_universe_Inflation_rate = Inflation_rates[i] + Best_universe = numpy.array(Universes[i, :]) + + sorted_Inflation_rates = numpy.sort(Inflation_rates) + sorted_indexes = numpy.argsort(Inflation_rates) + + for newindex in range(0, N): + Sorted_universes[newindex, :] = numpy.array( + Universes[sorted_indexes[newindex], :] + ) + + normalized_sorted_Inflation_rates = numpy.copy(normr(sorted_Inflation_rates)) + + Universes[0, :] = numpy.array(Sorted_universes[0, :]) + + for i in range(1, N): + Back_hole_index = i + for j in range(0, dim): + r1 = random.random() + + if r1 < normalized_sorted_Inflation_rates[i]: + White_hole_index = RouletteWheelSelection(-sorted_Inflation_rates) + + if White_hole_index == -1: + White_hole_index = 0 + White_hole_index = 0 + Universes[Back_hole_index, j] = Sorted_universes[ + White_hole_index, j + ] + + r2 = random.random() + + if r2 < WEP: + r3 = random.random() + if r3 < 0.5: + Universes[i, j] = Best_universe[j] + TDR * ( + (ub[j] - lb[j]) * random.random() + lb[j] + ) # random.uniform(0,1)+lb); + if r3 > 0.5: + Universes[i, j] = Best_universe[j] - TDR * ( + (ub[j] - lb[j]) * random.random() + lb[j] + ) # random.uniform(0,1)+lb); + + convergence[Time - 1] = Best_universe_Inflation_rate + if Time % 1 == 0: + print( + [ + "At iteration " + + str(Time) + + " the best fitness is " + + str(Best_universe_Inflation_rate) + ] + ) + + Time = Time + 1 + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence + s.optimizer = "MVO" + s.bestIndividual = Best_universe + s.objfname = objf.__name__ + + return s diff --git a/evolopy/PSO.py b/evolopy/PSO.py new file mode 100644 index 0000000..34876bf --- /dev/null +++ b/evolopy/PSO.py @@ -0,0 +1,103 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun May 15 22:37:00 2016 + +@author: Hossam Faris +""" + +import random +import numpy +from solution import solution +import time + + +def PSO(objf, lb, ub, dim, PopSize, iters): + + # PSO parameters + + Vmax = 6 + wMax = 0.9 + wMin = 0.2 + c1 = 2 + c2 = 2 + + s = solution() + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + ######################## Initializations + + vel = numpy.zeros((PopSize, dim)) + + pBestScore = numpy.zeros(PopSize) + pBestScore.fill(float("inf")) + + pBest = numpy.zeros((PopSize, dim)) + gBest = numpy.zeros(dim) + + gBestScore = float("inf") + + pos = numpy.zeros((PopSize, dim)) + for i in range(dim): + pos[:, i] = numpy.random.uniform(0, 1, PopSize) * (ub[i] - lb[i]) + lb[i] + + convergence_curve = numpy.zeros(iters) + + ############################################ + print('PSO is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + for l in range(0, iters): + for i in range(0, PopSize): + # pos[i,:]=checkBounds(pos[i,:],lb,ub) + for j in range(dim): + pos[i, j] = numpy.clip(pos[i, j], lb[j], ub[j]) + # Calculate objective function for each particle + fitness = objf(pos[i, :]) + + if pBestScore[i] > fitness: + pBestScore[i] = fitness + pBest[i, :] = pos[i, :].copy() + + if gBestScore > fitness: + gBestScore = fitness + gBest = pos[i, :].copy() + + # Update the W of PSO + w = wMax - l * ((wMax - wMin) / iters) + + for i in range(0, PopSize): + for j in range(0, dim): + r1 = random.random() + r2 = random.random() + vel[i, j] = ( + w * vel[i, j] + + c1 * r1 * (pBest[i, j] - pos[i, j]) + + c2 * r2 * (gBest[j] - pos[i, j]) + ) + + if vel[i, j] > Vmax: + vel[i, j] = Vmax + + if vel[i, j] < -Vmax: + vel[i, j] = -Vmax + + pos[i, j] = pos[i, j] + vel[i, j] + + convergence_curve[l] = gBestScore + + if l % 1 == 0: + print(["At iteration " + str(l + 1) + " the best fitness is " + str(gBestScore)]) + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence_curve + s.optimizer = "PSO" + s.bestIndividual = gBest + s.objfname = objf.__name__ + + return s diff --git a/evolopy/SCA.py b/evolopy/SCA.py new file mode 100644 index 0000000..2011c74 --- /dev/null +++ b/evolopy/SCA.py @@ -0,0 +1,93 @@ +""" Sine Cosine OPtimization Algorithm """ + +import random +import numpy +import math +from solution import solution +import time + + +def SCA(objf, lb, ub, dim, SearchAgents_no, Max_iter): + + # destination_pos + Dest_pos = numpy.zeros(dim) + Dest_score = float("inf") + + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # Initialize the positions of search agents + Positions = numpy.zeros((SearchAgents_no, dim)) + for i in range(dim): + Positions[:, i] = ( + numpy.random.uniform(0, 1, SearchAgents_no) * (ub[i] - lb[i]) + lb[i] + ) + + Convergence_curve = numpy.zeros(Max_iter) + s = solution() + + # Loop counter + print('SCA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + # Main loop + for l in range(0, Max_iter): + for i in range(0, SearchAgents_no): + + # Return back the search agents that go beyond the boundaries of the search space + for j in range(dim): + Positions[i, j] = numpy.clip(Positions[i, j], lb[j], ub[j]) + + # Calculate objective function for each search agent + fitness = objf(Positions[i, :]) + + if fitness < Dest_score: + Dest_score = fitness # Update Dest_Score + Dest_pos = Positions[i, :].copy() + + # Eq. (3.4) + a = 2 + Max_iteration = Max_iter + r1 = a - l * ((a) / Max_iteration) # r1 decreases linearly from a to 0 + + # Update the Position of search agents + for i in range(0, SearchAgents_no): + for j in range(0, dim): + + # Update r2, r3, and r4 for Eq. (3.3) + r2 = (2 * numpy.pi) * random.random() + r3 = 2 * random.random() + r4 = random.random() + + # Eq. (3.3) + if r4 < (0.5): + # Eq. (3.1) + Positions[i, j] = Positions[i, j] + ( + r1 * numpy.sin(r2) * abs(r3 * Dest_pos[j] - Positions[i, j]) + ) + else: + # Eq. (3.2) + Positions[i, j] = Positions[i, j] + ( + r1 * numpy.cos(r2) * abs(r3 * Dest_pos[j] - Positions[i, j]) + ) + + Convergence_curve[l] = Dest_score + + if l % 1 == 0: + print( + ["At iteration " + str(l) + " the best fitness is " + str(Dest_score)] + ) + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "SCA" + s.bestIndividual = Dest_pos + s.objfname = objf.__name__ + + return s diff --git a/evolopy/SSA.py b/evolopy/SSA.py new file mode 100644 index 0000000..8b880fe --- /dev/null +++ b/evolopy/SSA.py @@ -0,0 +1,125 @@ +import random +import numpy +import math +from solution import solution +import time + + +def SSA(objf, lb, ub, dim, N, Max_iteration): + + # Max_iteration=1000 + # lb=-100 + # ub=100 + # dim=30 + N = 50 # Number of search agents + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + Convergence_curve = numpy.zeros(Max_iteration) + + # Initialize the positions of salps + SalpPositions = numpy.zeros((N, dim)) + for i in range(dim): + SalpPositions[:, i] = numpy.random.uniform(0, 1, N) * (ub[i] - lb[i]) + lb[i] + SalpFitness = numpy.full(N, float("inf")) + + FoodPosition = numpy.zeros(dim) + FoodFitness = float("inf") + # Moth_fitness=numpy.fell(float("inf")) + + s = solution() + + print('SSA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + + for i in range(0, N): + # evaluate moths + SalpFitness[i] = objf(SalpPositions[i, :]) + + sorted_salps_fitness = numpy.sort(SalpFitness) + I = numpy.argsort(SalpFitness) + + Sorted_salps = numpy.copy(SalpPositions[I, :]) + + FoodPosition = numpy.copy(Sorted_salps[0, :]) + FoodFitness = sorted_salps_fitness[0] + + Iteration = 1 + + # Main loop + while Iteration < Max_iteration: + + # Number of flames Eq. (3.14) in the paper + # Flame_no=round(N-Iteration*((N-1)/Max_iteration)); + + c1 = 2 * math.exp(-((4 * Iteration / Max_iteration) ** 2)) + # Eq. (3.2) in the paper + + for i in range(0, N): + + SalpPositions = numpy.transpose(SalpPositions) + + if i < N / 2: + for j in range(0, dim): + c2 = random.random() + c3 = random.random() + # Eq. (3.1) in the paper + if c3 < 0.5: + SalpPositions[j, i] = FoodPosition[j] + c1 * ( + (ub[j] - lb[j]) * c2 + lb[j] + ) + else: + SalpPositions[j, i] = FoodPosition[j] - c1 * ( + (ub[j] - lb[j]) * c2 + lb[j] + ) + + #################### + + elif i >= N / 2 and i < N + 1: + point1 = SalpPositions[:, i - 1] + point2 = SalpPositions[:, i] + + SalpPositions[:, i] = (point2 + point1) / 2 + # Eq. (3.4) in the paper + + SalpPositions = numpy.transpose(SalpPositions) + + for i in range(0, N): + + # Check if salps go out of the search spaceand bring it back + for j in range(dim): + SalpPositions[i, j] = numpy.clip(SalpPositions[i, j], lb[j], ub[j]) + + SalpFitness[i] = objf(SalpPositions[i, :]) + + if SalpFitness[i] < FoodFitness: + FoodPosition = numpy.copy(SalpPositions[i, :]) + FoodFitness = SalpFitness[i] + + # Display best fitness along the iteration + if Iteration % 1 == 0: + print( + [ + "At iteration " + + str(Iteration) + + " the best fitness is " + + str(FoodFitness) + ] + ) + + Convergence_curve[Iteration] = FoodFitness + + Iteration = Iteration + 1 + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = Convergence_curve + s.optimizer = "SSA" + s.bestIndividual = FoodPosition + s.objfname = objf.__name__ + + return s diff --git a/evolopy/WOA.py b/evolopy/WOA.py new file mode 100644 index 0000000..8f784ef --- /dev/null +++ b/evolopy/WOA.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon May 16 14:19:49 2016 + +@author: hossam +""" +import random +import numpy +import math +from solution import solution +import time + + +def WOA(objf, lb, ub, dim, SearchAgents_no, Max_iter): + + # dim=30 + # SearchAgents_no=50 + # lb=-100 + # ub=100 + # Max_iter=500 + if not isinstance(lb, list): + lb = [lb] * dim + if not isinstance(ub, list): + ub = [ub] * dim + + # initialize position vector and score for the leader + Leader_pos = numpy.zeros(dim) + Leader_score = float("inf") # change this to -inf for maximization problems + + # Initialize the positions of search agents + Positions = numpy.zeros((SearchAgents_no, dim)) + for i in range(dim): + Positions[:, i] = ( + numpy.random.uniform(0, 1, SearchAgents_no) * (ub[i] - lb[i]) + lb[i] + ) + + # Initialize convergence + convergence_curve = numpy.zeros(Max_iter) + + ############################ + s = solution() + + print('WOA is optimizing "' + objf.__name__ + '"') + + timerStart = time.time() + s.startTime = time.strftime("%Y-%m-%d-%H-%M-%S") + ############################ + + t = 0 # Loop counter + + # Main loop + while t < Max_iter: + for i in range(0, SearchAgents_no): + + # Return back the search agents that go beyond the boundaries of the search space + + # Positions[i,:]=checkBounds(Positions[i,:],lb,ub) + for j in range(dim): + Positions[i, j] = numpy.clip(Positions[i, j], lb[j], ub[j]) + + # Calculate objective function for each search agent + fitness = objf(Positions[i, :]) + + # Update the leader + if fitness < Leader_score: # Change this to > for maximization problem + Leader_score = fitness + # Update alpha + Leader_pos = Positions[ + i, : + ].copy() # copy current whale position into the leader position + + a = 2 - t * ((2) / Max_iter) + # a decreases linearly fron 2 to 0 in Eq. (2.3) + + # a2 linearly decreases from -1 to -2 to calculate t in Eq. (3.12) + a2 = -1 + t * ((-1) / Max_iter) + + # Update the Position of search agents + for i in range(0, SearchAgents_no): + r1 = random.random() # r1 is a random number in [0,1] + r2 = random.random() # r2 is a random number in [0,1] + + A = 2 * a * r1 - a # Eq. (2.3) in the paper + C = 2 * r2 # Eq. (2.4) in the paper + + b = 1 + # parameters in Eq. (2.5) + l = (a2 - 1) * random.random() + 1 # parameters in Eq. (2.5) + + p = random.random() # p in Eq. (2.6) + + for j in range(0, dim): + + if p < 0.5: + if abs(A) >= 1: + rand_leader_index = math.floor( + SearchAgents_no * random.random() + ) + X_rand = Positions[rand_leader_index, :] + D_X_rand = abs(C * X_rand[j] - Positions[i, j]) + Positions[i, j] = X_rand[j] - A * D_X_rand + + elif abs(A) < 1: + D_Leader = abs(C * Leader_pos[j] - Positions[i, j]) + Positions[i, j] = Leader_pos[j] - A * D_Leader + + elif p >= 0.5: + + distance2Leader = abs(Leader_pos[j] - Positions[i, j]) + # Eq. (2.5) + Positions[i, j] = ( + distance2Leader * math.exp(b * l) * math.cos(l * 2 * math.pi) + + Leader_pos[j] + ) + + convergence_curve[t] = Leader_score + if t % 1 == 0: + print( + ["At iteration " + str(t) + " the best fitness is " + str(Leader_score)] + ) + t = t + 1 + + timerEnd = time.time() + s.endTime = time.strftime("%Y-%m-%d-%H-%M-%S") + s.executionTime = timerEnd - timerStart + s.convergence = convergence_curve + s.optimizer = "WOA" + s.objfname = objf.__name__ + s.best = Leader_score + s.bestIndividual = Leader_pos + + return s