Genetic Algorithm Traveling Salesman Problem Solver

tsp_ga.py

import math
import random
import numpy as np
from multiprocessing import Pool, cpu_count

def prepend_start(route):
    return [0] + list(route)


def create_population(num_points, size):
    order = list(range(1, num_points))
    
    return [np.random.permutation(order) for _ in range(size)]


def get_pool(n_jobs):
    num_processes = int(n_jobs if n_jobs > 0 else cpu_count())
    
    if num_processes > 1:
        return Pool(num_processes)
        

def compute_fitness(fitness_fn, points, routes, pool=None):    
    if pool is None:
        return [fitness_fn(points, prepend_start(route)) for route in routes]    
    
    return pool.starmap(fitness_fn, ((points, prepend_start(route)) for route in routes))


def select_parents(fitness):
    indices = np.argsort(fitness)
    
    return indices[:int(len(fitness)*0.5)]
    

def crossover(routes, parents, target_size):
    random.shuffle(parents)
    
    num_children = 2 * int(math.ceil(1.0*target_size / len(parents)))
    
    new_routes = []
    for i in range(0, len(parents)-1, 2):
        parentA, parentB = routes[parents[i]], routes[parents[i+1]]
        
        for _ in range(num_children):
            child = np.empty(len(parentA), dtype=np.int8)
            
            split = np.random.randint(0, int(len(parentA) / 2.0))
            span = (split, split + int(len(parentA) / 4.0) + 1)
            
            child[span[0]:span[1]] = parentA[span[0]:span[1]]
            
            free_indices = list(range(0, span[0])) + list(range(span[1], len(parentA)))
            used_values = set(parentA[span[0]:span[1]])
            
            i = 0
            for point in parentB:
                if point not in used_values:
                    child[free_indices[i]] = point
                    i += 1
                    
            new_routes.append(child)
    
    return new_routes    


def mutate(routes, p=0.1):
    for i, route in enumerate(routes):
        if np.random.uniform() <= p:
            idxA = np.random.randint(len(route))
            idxB = idxA
            while idxB == idxA:
                idxB = np.random.randint(len(route))
                
            route[idxA], route[idxB] = route[idxB], route[idxA]
            routes[i] = route
            
    return routes


def optimize_route(points, fitness_fn, min_population_size=100, max_population_size=1000, max_iter=1000, mutation_p=0.1, n_jobs=-1):
    num_points = len(points)
    population_size = int(np.clip(num_points * 10, min_population_size, max_population_size))
    
    pool = get_pool(n_jobs)
    
    population = create_population(num_points, population_size)
    fitness_means = []
    for i in range(max_iter):
        fitness = compute_fitness(fitness_fn, points, population, pool)
        parents = select_parents(fitness)
        population = crossover(population, parents, population_size)
        population = mutate(population, p=mutation_p)
        
        fitness_means.append(np.mean(fitness))
        
        if i > 0 and (i % 50) == 0:
            threshold = fitness_means[-1] * 0.001
            if np.abs(np.mean(fitness_means[-50:]) - np.mean(fitness_means[-10:])) < threshold:
                break
        
    fitness = compute_fitness(fitness_fn, points, population, pool)
    best_idx = np.argmin(fitness)
    
    if pool is not None:
        pool.close()
        pool.join()
    
    return prepend_start(population[best_idx])