# Genetic algorithm to guess coefficient of a polynomial

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
+ array * x
+ array * x ** 2
+ array * x ** 3
+ array * 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
+ ind * x
+ ind * x ** 2
+ ind * x ** 3
+ ind * x ** 4
)
- y
) ** 2
ind = 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] < point < intervals[j]:
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)
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])

$$$$


# This is not the best algorithm

If the goal is to get the best coefficients for a polynomial so it fits the given points, then a polynomial regression algorithm such as numpy.polynomial.polynomial.Polynomial.fit() will give you the best fit much faster, as there is an analytic solution to the polynomial least squares problem.

If the goal is to learn about genetic algorithms, you can ignore this though.

# Avoid hardcoding numbers

Hardcoding numbers should be avoided where possible. For example, the initial array in the fitness() function can be derived from poly():

initial = [[i, poly(i)] for i in range(-5, 6)]


Some parts of your code are unnecessarily complicated. For example, when creating prob_list in sel_ind(), you create a deep copy of the input array, then only modify parts of it, then create a list from the modified parts. This can be simplified a lot by first just taking the elements you want, then normalizing them, without needing to make a deep copy or force the creation of a list:

ranks = nk_array[:, -1]
norm = sum(ranks)
probabilities = ranks / norm
...
for prob in probabilities:
...


Generating the list of intervals can also be simplified using this trick with zip():

intervals = list(zip( + probabilities, probabilities))


List comprehension can also simplify some code (although I recommend you only use it for simple expressions, as list comprehension can quickly get... incomprehensible):

start = rng.uniform(0, 1 / 5)
points = [start + i / 5 for i in range(5)]


# Avoid code duplication

Related to the above, whenever you find you are repeating things, it's time to create a loop or a function to avoid code duplication. In crossover(), you are creating two "babies", but instead of writing out the code for each baby separately, find a way to write the code only once. I would just create a function here (it can be nested inside the for-loops):

def make_baby(array, i, j):
mate = rng.uniform(-0.25, 1.25, [1,5]).squeeze()
mutation = rng.uniform(-0.025, 0.025, [1, 5]).squeeze()
return np.append((array[i, :5] * mate + array[j, :5] * (1 - mate) + mutation * mate), 0)

progeny = np.append(progeny, [make_baby(array, i, j)], axis=0)
progeny = np.append(progeny, [make_baby(array, j, i)], axis=0)


In this case you could also have avoided the need for duplication by having i and j both loop over the full range(len(individuals)), and then just have an if i == j: continue inside the loop.

# Split the code into even more functions

The __main__ part is still quite large, and could be split up to make it more concise and readable. Ideally, it should be a short function that gives you a high level overview of what your program is doing. It could look like this:

if __name__ == "__main__":
# Initialize
rng = np.random.default_rng()
fig, ax = initialize_plot()
pool = initialize_pool()
print("Starting fitness: ", pool[0, 5])

# Run the loop until we converge
loop_count = 0

while pool[0, 5] > 0.005:
if loop_count % 100 == 0:
update_plot(fig, ax, pool)

update_pool(pool)
loop_count += 1

# Show the final result
update_plot(fig, ax, pool)
print("Final fitness: ", pool[0, 5])


# Naming things

Try to be a bit more precise when naming things, and avoid unnecessary abbreviations. For example, sel_ind() could be renamed to select_individuals(), which would make it much clearer what this function does, without having to read the comments above the function definition.

Also avoid reusing the same variable name for different things, for example like you did in these two lines:

pool = rng.uniform(-50, 50, [50, 6])
pool = fitness(pool)


Either avoid the reuse, like so:

pool = fitness(rng.uniform(-50, 50, [50, 6]))


Or give the first use a better name:

random_coefficients = rng.uniform(...)
pool = fitness(random_coefficients)


I would also use verbs for functions, and nouns for variables and constants. For variables, use the plural if it is a colection of values, like a list or array. There are exceptions to these rules, for example if the word itself already suggests it is a function (like poly()) or a collection of things (for an obvious example, something named array). So I suggest:

• initial -> initial_points
• fitness() -> calculate_fitness()
• nl_rank() -> calculate_nl_rank()
• sel_ind() -> select_individuals()
• index -> indices (in the above function)
• crossover() -> generate_crossovers()
• arr_print() -> print_array()
• My goal was to learn about genetic algorithms and not using the most optimal method available. I will try to make changes according to your suggestions. – Roni Saiba Jun 12 at 9:56
• I do not know enough about genetic algorithms, but are we sure that the final while loop converges? There is a terminate condition on elif pool[0, 5] < 0.005: but no terminate condition on count_loops`. – N3buchadnezzar Jun 12 at 10:41
• @N3buchadnezzar For the conditions I used it does converge. I first ran it with termination on loop count to get an idea of how fast convergence is. Then I tried pool [0,5] < 10, < 5 and successively smaller numebrs. – Roni Saiba Jun 12 at 13:08