Genetic algorithm to guess coefficient of a polynomial

Question

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])

```

Answer 1

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)]

Simplify your code

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([0] + 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()