I have tried to code a genetic algorithm to guess the coefficients of a degree 4 polynomial. The information initially provided is values of y = f(x) for different x using the original polynomial. I then generate 100 polynomials with randomly selected coefficients. These polynomials are then ranked based on the least square difference (LSD) and sorted such that lower LSD has higher rank. These ranks are then converted into nonlinear ranking and mapped to a line of unit length. Individuals are then selected and real valued crossover is used to generate progeny, who are also sorted according to LSD. The worst individuals in the starting pool are then replaced by progeny. Please suggest improvements where possible. I am intermediate in Python.
"""
Genetic Algorithm implementation of finding coefficients of a polynomial
3.0 - 4.3 * x + 5.9 * x ** 2 - 5.2 * x ** 3 + 1.0 * x ** 4 = 0
"""
#%% import modules
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
# np.set_printoptions(linewidth=np.inf)
#%% define a function to return the polynomial
def poly(x):
return 3.0 - 4.3 * x + 5.9 * x ** 2 - 5.2 * x ** 3 + 1.0 * x ** 4
def gen_poly(array, x):
return (
array[0]
+ array[1] * x
+ array[2] * x ** 2
+ array[3] * x ** 3
+ array[4] * x ** 4
)
#%% define function to measure fitness
### takes an array as input and fills the last column with fitness values
### take care to always keep one column for fitness values in starting array
### output is an array sorted by fitness
def fitness(array):
# set initial values
initial = np.array(
[
[-5, 1447.0],
[-4, 703.4],
[-3, 290.4],
[-2, 92.8],
[-1, 19.4],
[0, 3],
[1, 0.4],
[2, -7.6],
[3, -16.2],
[4, 3.4],
[5, 104],
]
)
# score is the total square deviation from initial. lower is better
for ind in array:
fitness = 0
for x, y in zip(initial[:, 0], initial[:, 1]):
fitness += (
(
ind[0]
+ ind[1] * x
+ ind[2] * x ** 2
+ ind[3] * x ** 3
+ ind[4] * x ** 4
)
- y
) ** 2
ind[5] = fitness
return array[array[:, 5].argsort()]
#%% function to give non linear rank to sorted pool
def nl_rank(array):
copy = deepcopy(array)
# normalization factor for selection pressure 3 and individuals 50
norm = 290.3359045831912
# setting ranks based on position in the pool
for i, j in enumerate(copy):
j[-1] = 50 * 1.06 ** (49 - i) / norm
return copy
#%% define function to carry out stochastic universal sampling
### takes non linear ranked array as input
### output is list of selected individuals
def sel_ind(nl_array):
copy = deepcopy(nl_array)
# normalize ranks
norm = sum(copy[:, -1])
copy[:, -1] = copy[:, -1] / norm
# map intervals on range [0, 1]
prob_list = list(copy[:, -1])
intervals = []
start = 0
for prob in prob_list:
end = start + prob
intervals.append((start, end))
start = end
# selecting 6 individuals from the intervals
rng = np.random.default_rng()
points = [rng.uniform(0, 1 / 5)]
for i in range(4):
points.append(points[-1] + 1 / 5)
index, i = [], 0
for point in points:
for j in range(i, len(intervals)):
if intervals[j][0] < point < intervals[j][1]:
index.append(j)
i = j
break
return index
#%% define function to carry out mating. only unique pairings are considered
### each mating gives 2 children
def crossover(array, individuals):
rng = np.random.default_rng()
progeny = np.empty((0, 6))
for i in range(len(individuals) - 1):
for j in range(i + 1, len(individuals)):
mate1 = rng.uniform(-0.25, 1.25, [1, 5]).squeeze()
mate2 = rng.uniform(-0.25, 1.25, [1, 5]).squeeze()
mutation1 = rng.uniform(-0.025, 0.025, [1, 5]).squeeze()
mutation2 = rng.uniform(-0.025, 0.025, [1, 5]).squeeze()
baby1 = (
array[i, :5] * mate1 + array[j, :5] * (1 - mate1) + mutation1 * mate1
)
baby1 = np.append(baby1, 0)
progeny = np.append(progeny, [baby1], axis=0)
baby2 = (
array[j, :5] * mate2 + array[i, :5] * (1 - mate2) + mutation2 * mate2
)
baby2 = np.append(baby2, 0)
progeny = np.append(progeny, [baby2], axis=0)
return fitness(progeny)
#%% helper function to print arrays to log
def arr_print(arr, count):
print(f"#loop_count = {count}")
print(arr)
#%% main
if __name__ == "__main__":
# create rng instance
rng = np.random.default_rng()
# create parent pool
# parent pool has 5 columns for coefficients and one for fitness
pool = rng.uniform(-50, 50, [50, 6])
# measure fitness of each parent and sort in decreasing order of fitness
pool = fitness(pool)
starting_fitness = pool[0, 5]
# plotting the original curve
plt.ion()
fig, ax = plt.subplots()
x = np.linspace(-5, 5, 200)
ax.plot(x, poly(x), "r", lw=1, label="Original")
ax.set_xlabel("X axis")
ax.set_ylabel("Y axis")
ax.set_ylim(-1000, 1000)
ax.set_title("Comparision of curves")
ax.legend()
loop_count = 1
while True:
# plotting
if loop_count == 1:
ax.plot(x, gen_poly(pool[0, :], x), lw=1, label=f"Iteration = {loop_count}")
ax.legend()
fig.canvas.draw()
plt.pause(0.0001)
# arr_print(pool, loop_count)
elif loop_count % 100 == 0:
ax.plot(x, gen_poly(pool[0, :], x), lw=0.25, ls="solid")
fig.canvas.draw()
plt.pause(0.0000001)
# arr_print(pool, loop_count)
# add termination condition
elif pool[0, 5] < 0.005:
ax.plot(
x,
gen_poly(pool[0, :], x),
"k",
lw=1.5,
label=f"Iteration = {loop_count}",
)
ax.set_xlim(-2, 5)
ax.set_ylim(-40, 75)
ax.legend()
fig.canvas.draw()
plt.pause(0.001)
# arr_print(pool, loop_count)
break
# rank parents based on non linear ranking
ranked = nl_rank(pool)
# select individuals
individuals = sel_ind(ranked)
# create progeny
progeny = crossover(ranked, individuals)
# remove 20 worst individuals from pool
pool = np.delete(pool, np.s_[-20:], axis=0)
# add progeny to the new pool
pool = np.vstack((pool, progeny[:20, :]))
# sort pool according to fitness
pool = pool[pool[:, 5].argsort()]
loop_count += 1
print(starting_fitness, pool[0, 5])
```