Genetic Algorithms › Implementation in NumPy

import time
from itertools import permutations

import numpy as np
import pandas as pd
from book_plot_style import seed_book_plots

seed_book_plots()

BUDGET = 300
SEEDS = list(range(20))
CUSTOMERS = tuple(range(1, 9))

COORDS = np.array(
    [
        [0.0, 0.0],  # depot
        [1.0, 3.0],  # 1
        [2.0, 6.0],  # 2
        [4.0, 5.0],  # 3
        [7.0, 6.0],  # 4
        [8.0, 3.0],  # 5
        [6.0, 1.0],  # 6
        [3.0, 1.0],  # 7
        [5.0, 3.0],  # 8
    ]
)
N_NODES = len(COORDS)
DEPOT = 0
DEMAND = np.array([0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
CAPACITY = 4.0
WINDOW_OPEN = np.zeros(N_NODES)
WINDOW_CLOSE = np.full(N_NODES, 80.0)
SERVICE = np.array([0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

_pairs = COORDS[:, None, :] - COORDS[None, :, :]
DIST = np.sqrt((_pairs**2).sum(axis=2))

def decode(perm):
    """Split an 8-customer permutation into two four-customer depot tours."""
    values = np.asarray(perm)
    if (
        values.shape != (8,)
        or values.dtype.kind not in "iu"
        or sorted(values.tolist()) != list(CUSTOMERS)
    ):
        raise ValueError("Expected one integer permutation of customers1..8")
    perm = values.tolist()
    return [[DEPOT, *perm[0:4], DEPOT], [DEPOT, *perm[4:8], DEPOT]]


def schedule(route):
    """Simulate one depot tour: arc lengths, load, service starts, window checks."""
    load = 0.0
    clock = 0.0
    length = 0.0
    feasible = True
    for a, b in zip(route[:-1], route[1:]):
        load += float(DEMAND[a])
        travel = float(DIST[a, b])
        length += travel
        arrival = clock + float(SERVICE[a]) + travel
        start = max(arrival, float(WINDOW_OPEN[b]))
        if start > float(WINDOW_CLOSE[b]) + 1e-9:
            feasible = False
        clock = start
    if load > CAPACITY + 1e-9:
        feasible = False
    return length, load, clock, feasible


def objective(perm):
    """Complete-route objective: summed Euclidean length of both closed tours."""
    total = 0.0
    for route in decode(perm):
        length, _, _, feasible = schedule(route)
        if not feasible:
            return float("inf")
        # Match the preceding chapter's sequential10-arc accumulation.
        for a, b in zip(route[:-1], route[1:]):
            total += float(DIST[a, b])
    return total


def validate_solution(perm):
    """Independent check that recomputes distances from coordinates, not DIST."""
    values = np.asarray(perm)
    if (
        values.shape != (8,)
        or values.dtype.kind not in "iu"
        or sorted(values.tolist()) != list(CUSTOMERS)
    ):
        return {
            "routes": [],
            "total": float("inf"),
            "feasible": False,
            "permutation": values.tolist(),
        }
    perm = values.tolist()
    # Reconstruct the two routes independently of decode and DIST.
    routes = [[0, *perm[:4], 0], [0, *perm[4:], 0]]
    rows = []
    total = 0.0
    feasible = sorted(perm) == list(CUSTOMERS)
    for route in routes:
        load = 0.0
        clock = 0.0
        length = 0.0
        windows_ok = True
        for a, b in zip(route[:-1], route[1:]):
            load += float(DEMAND[a])
            travel = float(
                np.hypot(
                    COORDS[a, 0] - COORDS[b, 0], COORDS[a, 1] - COORDS[b, 1]
                )
            )
            length += travel
            arrival = clock + float(SERVICE[a]) + travel
            start = max(arrival, float(WINDOW_OPEN[b]))
            if start > float(WINDOW_CLOSE[b]) + 1e-9:
                windows_ok = False
            clock = start
        capacity_ok = load <= CAPACITY + 1e-9
        rows.append(
            {
                "route": route,
                "load": load,
                "length": length,
                "return_to_depot": clock,
                "capacity_ok": capacity_ok,
                "windows_ok": windows_ok,
            }
        )
        total += length
        feasible = feasible and capacity_ok and windows_ok
    return {
        "routes": rows,
        "total": total,
        "feasible": feasible,
        "permutation": perm,
    }

class BudgetedObjective:
    """Counts every complete-route evaluation and records the incumbent trace."""

    def __init__(self, budget):
        if (
            isinstance(budget, (bool, np.bool_))
            or not isinstance(budget, (int, np.integer))
            or budget < 1
        ):
            raise ValueError("budget must be a positive integer")
        self.budget = int(budget)
        self.used = 0
        self.best = float("inf")
        self.trace = []

    def __call__(self, perm):
        if self.used >= self.budget:
            raise RuntimeError("Objective evaluation budget exhausted")
        self.used += 1
        value = objective(perm)
        if value < self.best:
            self.best = value
        self.trace.append(self.best)
        return value

def nearest_neighbor_perm():
    """Greedy construction: nearest remaining customer within each vehicle block."""
    perm = []
    remaining = list(CUSTOMERS)
    current = DEPOT
    while remaining:
        if len(perm) % 4 == 0:
            current = DEPOT
        nxt = min(remaining, key=lambda c: (float(DIST[current, c]), c))
        perm.append(nxt)
        remaining.remove(nxt)
        current = nxt
    return np.array(perm, dtype=int)


NN_PERM = nearest_neighbor_perm()
NN_COST = float(objective(NN_PERM))

def nn_seeded_random_search(rng, budget, nn_perm):
    """Nearest-neighbour score, then uniform draws to the remaining budget."""
    ev = BudgetedObjective(budget)
    initial_cost = ev(nn_perm)

    best_perm = np.asarray(nn_perm, dtype=int).copy()
    best_cost = float(initial_cost)
    while ev.used < budget:
        perm = rng.permutation(CUSTOMERS)
        value = ev(perm)
        if value < best_cost:
            best_perm, best_cost = perm.copy(), value
    return best_perm, best_cost, ev


def hill_climbing(rng, budget, nn_perm):
    """Tolerance-aware scans of all 28 segment reversals, with restarts."""
    ev = BudgetedObjective(budget)
    initial_cost = ev(nn_perm)

    current = np.asarray(nn_perm, dtype=int).copy()
    current_cost = float(initial_cost)
    best_perm, best_cost = current.copy(), current_cost

    while ev.used < budget:
        neighbor_cost = float("inf")
        neighbor = None
        completed = True
        for i in range(8):
            for j in range(i + 1, 8):
                if ev.used >= budget:
                    completed = False
                    break
                trial = current.copy()
                trial[i : j + 1] = trial[i : j + 1][::-1]
                value = ev(trial)
                if (
                    value < neighbor_cost - 1e-12
                ):  # retain first on numerical ties
                    neighbor_cost = value
                    neighbor = trial
                if value < best_cost:
                    best_perm, best_cost = trial.copy(), value
            if not completed:
                break
        if not completed:
            break  # partial sweep at exhaustion
        if neighbor_cost < current_cost - 1e-12:
            current, current_cost = neighbor, neighbor_cost
        else:
            if ev.used >= budget:
                break
            current = rng.permutation(CUSTOMERS)  # restart, one counted call
            current_cost = ev(current)
            if current_cost < best_cost:
                best_perm, best_cost = current.copy(), current_cost
    return best_perm, best_cost, ev


def simulated_annealing(rng, budget, nn_perm, t_start=2.0, t_end=0.01):
    """Symmetric single-reversal proposal, geometric cooling, best-so-far reported."""
    if not (
        np.isfinite(t_start) and np.isfinite(t_end) and 0 < t_end <= t_start
    ):
        raise ValueError("Require finite temperatures with0<t_end<=t_start")
    temperature_ratio = t_end / t_start
    if temperature_ratio < np.finfo(np.asarray(temperature_ratio).dtype).tiny:
        raise ValueError(
            "temperature ratio must be normal in its floating-point dtype"
        )
    ev = BudgetedObjective(budget)
    initial_cost = ev(nn_perm)

    current = np.asarray(nn_perm, dtype=int).copy()
    current_cost = float(initial_cost)
    best_perm, best_cost = current.copy(), current_cost

    while ev.used < budget:
        fraction = ev.used / max(budget - 1, 1)
        temperature = t_start * (t_end / t_start) ** fraction
        i, j = sorted(int(v) for v in rng.choice(8, size=2, replace=False))
        trial = current.copy()
        trial[i : j + 1] = trial[i : j + 1][::-1]
        value = ev(trial)
        delta = value - current_cost
        if delta < 0.0 or rng.random() < np.exp(-delta / temperature):
            current, current_cost = trial, value
            if value < best_cost:
                best_perm, best_cost = trial.copy(), value
    return best_perm, best_cost, ev

def order_crossover(parent_one, parent_two, a, b):
    """Order crossover: copy one block, fill the rest in cyclic order from parent two."""
    n = len(parent_one)
    child = np.full(n, -1, dtype=int)
    child[a : b + 1] = parent_one[a : b + 1]
    used = set(int(v) for v in child[a : b + 1])

    empty_positions = [(b + 1 + k) % n for k in range(n - (b - a + 1))]
    scan_order = [(b + 1 + k) % n for k in range(n)]
    values = [
        int(parent_two[t]) for t in scan_order if int(parent_two[t]) not in used
    ]

    for position, value in zip(empty_positions, values):
        child[position] = value
    return child


def swap_mutation(perm, rng):
    """Exchange the values at two distinct positions; validity is automatic."""
    child = perm.copy()
    i, j = (int(v) for v in rng.choice(len(child), size=2, replace=False))
    child[i], child[j] = child[j], child[i]
    return child


def tournament(rng, population, size):
    """Sample distinct individuals and return the permutation with lowest cost."""
    picked = rng.choice(len(population), size=size, replace=False)
    contenders = [population[int(k)] for k in picked]
    return min(contenders, key=lambda individual: individual[1])[0]


def elitist_survivors(parents, children, population_size):
    """Mu-plus-lambda survival: keep the cheapest `size` of parents plus children."""
    pool = parents + children
    pool = sorted(
        pool, key=lambda individual: individual[1]
    )  # stable: parents win ties
    return [(perm.copy(), float(cost)) for perm, cost in pool[:population_size]]


def canonical_signature(perm):
    """Vehicle-symmetric signature: tours are unordered and reversals are equal."""
    signature = []
    for route in decode(perm):
        sequence = tuple(int(v) for v in route[1:-1])
        signature.append(min(sequence, sequence[::-1]))
    return tuple(sorted(signature))


def population_stats(population):
    """Diversity and cost summary of one generation."""
    costs = np.array([cost for _, cost in population], dtype=float)
    matrix = np.array([perm for perm, _ in population], dtype=int)
    hamming = (matrix[:, None, :] != matrix[None, :, :]).sum(axis=2)
    upper = np.triu_indices(len(matrix), k=1)
    return {
        "best": float(costs.min()),
        "median": float(np.median(costs)),
        "worst": float(costs.max()),
        "distinct": int(
            len({canonical_signature(perm) for perm, _ in population})
        ),
        "mean_hamming": float(hamming[upper].mean()),
    }

def genetic_algorithm(
    rng,
    budget,
    nn_perm,
    pop_size=20,
    tournament_size=3,
    crossover_rate=0.9,
    mutation_rate=0.5,
):
    """Elitist generational GA over permutations. Returns (best_perm, best_cost, counter, history)."""
    n = len(CUSTOMERS)
    ev = BudgetedObjective(budget)
    if (
        isinstance(pop_size, (bool, np.bool_))
        or not isinstance(pop_size, (int, np.integer))
        or pop_size < 2
        or budget < pop_size
    ):
        raise ValueError(
            "Require population size>=2 and budget>=population size"
        )
    if (
        isinstance(tournament_size, (bool, np.bool_))
        or not isinstance(tournament_size, (int, np.integer))
        or not 1 <= tournament_size <= pop_size
    ):
        raise ValueError("Tournament size must be within the population")
    if not 0 <= crossover_rate <= 1 or not 0 <= mutation_rate <= 1:
        raise ValueError("Variation probabilities must lie in[0,1]")
    initial_cost = ev(nn_perm)
    population = [(np.asarray(nn_perm, dtype=int).copy(), initial_cost)]
    while len(population) < pop_size and ev.used < budget:
        perm = rng.permutation(CUSTOMERS)
        population.append((perm, ev(perm)))

    history = []
    while ev.used < budget:
        children = []
        for _ in range(pop_size):
            if ev.used >= budget:
                break
            parent_one = tournament(rng, population, tournament_size)
            parent_two = tournament(rng, population, tournament_size)
            if rng.random() < crossover_rate:
                a = int(rng.integers(0, n - 1))
                b = int(rng.integers(a + 1, n))
                child = order_crossover(parent_one, parent_two, a, b)
            else:
                child = parent_one.copy()
            if rng.random() < mutation_rate:
                child = swap_mutation(child, rng)
            children.append((child, ev(child)))
        population = elitist_survivors(population, children, pop_size)
        history.append(population_stats(population))

    best_perm, best_cost = min(population, key=lambda individual: individual[1])
    return best_perm.copy(), float(best_cost), ev, history

example_perms = {
    "P1 nearest neighbour": [1, 7, 8, 3, 6, 5, 4, 2],
    "P2": [1, 2, 3, 8, 6, 5, 4, 7],
    "P3": [1, 2, 3, 7, 4, 5, 6, 8],
    "P4": [2, 6, 1, 5, 3, 7, 4, 8],
}
child = order_crossover(
    np.array(example_perms["P2"]), np.array(example_perms["P3"]), a=4, b=6
)
mutated = np.array([1, 2, 3, 7, 6, 4, 5, 8])

example_table = pd.DataFrame(
    [
        {
            "individual": name,
            "permutation": tuple(perm),
            "cost": objective(perm),
        }
        for name, perm in example_perms.items()
    ]
    + [
        {
            "individual": "OX child",
            "permutation": tuple(int(v) for v in child),
            "cost": objective(child),
        },
        {
            "individual": "OX child + swap",
            "permutation": tuple(int(v) for v in mutated),
            "cost": objective(mutated),
        },
    ]
)

def run_ga(rng, budget, nn_perm):
    best_perm, best_cost, ev, _ = genetic_algorithm(rng, budget, nn_perm)
    return best_perm, best_cost, ev


RUNNERS = {
    "NN-seeded random search": nn_seeded_random_search,
    "2-opt hill climbing": hill_climbing,
    "Simulated annealing": simulated_annealing,
    "Genetic algorithm": run_ga,
}

records = []
traces = {name: [] for name in RUNNERS}
final_perms = {name: {} for name in RUNNERS}
ga_histories = {}

for seed in SEEDS:
    for name, runner in RUNNERS.items():
        rng = np.random.default_rng(seed)
        started = time.perf_counter()
        if name == "Genetic algorithm":
            best_perm, best_cost, ev, history = genetic_algorithm(
                rng, BUDGET, NN_PERM
            )
            ga_histories[seed] = history
        else:
            best_perm, best_cost, ev = runner(rng, BUDGET, NN_PERM)
        elapsed = time.perf_counter() - started
        records.append(
            {
                "method": name,
                "seed": seed,
                "best_cost": float(best_cost),
                "evaluations": int(ev.used),
                "seconds": elapsed,
            }
        )
        traces[name].append(np.asarray(ev.trace, dtype=float))
        final_perms[name][seed] = np.asarray(best_perm, dtype=int)

results = pd.DataFrame(records)
assert (results["evaluations"] == BUDGET).all(), "budget accounting violated"

all_perms = np.array(list(permutations(CUSTOMERS)), dtype=int)
exact_costs = np.array([objective(perm) for perm in all_perms])
EXACT_BEST = float(exact_costs.min())
exact_best_perms = all_perms[
    np.isclose(exact_costs, EXACT_BEST, rtol=0.0, atol=1e-10)
]
distinct_signatures = len({canonical_signature(perm) for perm in all_perms})

grouped = results.groupby("method")
summary = pd.DataFrame(
    {
        "mean": grouped["best_cost"].mean(),
        "std": grouped["best_cost"].std(ddof=1),
        "median": grouped["best_cost"].median(),
        "min": grouped["best_cost"].min(),
        "max": grouped["best_cost"].max(),
        "median_seconds": grouped["seconds"].median(),
    }
).reindex(list(RUNNERS))

summary["hits_exact_optimum"] = (
    results.assign(
        hit=np.isclose(results["best_cost"], EXACT_BEST, rtol=0.0, atol=1e-10)
    )
    .groupby("method")["hit"]
    .sum()
    .reindex(list(RUNNERS))
    .astype(int)
)
summary["mean_gap_pct"] = 100.0 * (summary["mean"] - EXACT_BEST) / EXACT_BEST

best_rows = results.loc[results.groupby("method")["best_cost"].idxmin()]
validation_rows = []
for _, row in results.iterrows():
    perm = final_perms[row["method"]][int(row["seed"])]
    report = validate_solution(perm)
    assert report["feasible"]
    assert np.isclose(report["total"], row["best_cost"], rtol=0.0, atol=1e-10)
    validation_rows.append(
        {
            "method": row["method"],
            "seed": int(row["seed"]),
            "tour 1": " - ".join(str(v) for v in report["routes"][0]["route"]),
            "tour 2": " - ".join(str(v) for v in report["routes"][1]["route"]),
            "length": report["total"],
            "feasible": report["feasible"],
            "matches report": bool(
                np.isclose(
                    report["total"], row["best_cost"], rtol=0.0, atol=1e-10
                )
            ),
        }
    )
validation_all = pd.DataFrame(validation_rows)
validation = validation_all.loc[best_rows.index].reset_index(drop=True)