Chapter 54: Combinatorial Optimization

import math
import random

import numpy as np

from kaggle_companion.activities._common import clean

SEED = 54
CITIES, BUDGET, INSTANCES = 150, 60000, 8   # every run gets the same 60,000 evaluated 2-opt moves
END_TEMPERATURE = 0.01                      # final temperature, as a fraction of the mean edge of the nearest-neighbour tour


def make_instance(seed):
    """Random Euclidean cities as a distance table (nested lists: fastest for scalar lookups)."""
    points = np.random.default_rng(seed).random((CITIES, 2))
    return np.sqrt(((points[:, None] - points[None]) ** 2).sum(-1)).tolist()


def tour_length(D, tour):
    return sum(D[tour[i]][tour[(i + 1) % CITIES]] for i in range(CITIES))


def nearest_neighbour_tour(D):
    """Cheap legal initial tour: always walk to the closest unvisited city."""
    left, tour = set(range(1, CITIES)), [0]
    while left:
        nxt = min(left, key=lambda j: D[tour[-1]][j])
        tour.append(nxt)
        left.remove(nxt)
    return tour


def two_opt_delta(D, tour, i, j):
    """Change in length if tour[i..j] is reversed (swap two edges)."""
    a, b, c, d = tour[i - 1], tour[i], tour[j], tour[(j + 1) % CITIES]
    return D[a][c] + D[b][d] - D[a][b] - D[c][d]


def sample_worsening(D, tour, rng, count=300):
    """Worsening 2-opt deltas sampled from the current tour."""
    out = []
    while len(out) < count:
        i, j = sorted(rng.sample(range(1, CITIES), 2))
        delta = two_opt_delta(D, tour, i, j)
        if delta > 0:
            out.append(delta)
    return out


def acceptance(deltas, temperature):
    """Mean probability that a sampled worsening move is accepted: mean of exp(-delta / T)."""
    return sum(math.exp(-d / temperature) for d in deltas) / len(deltas)


def solve_temperature(deltas, target):
    """The chapter's recipe: the acceptance is monotone in T, so bisect for the T that gives the target."""
    lo, hi = 1e-9, 1e3
    for _ in range(60):
        mid = math.sqrt(lo * hi)
        lo, hi = (mid, hi) if acceptance(deltas, mid) < target else (lo, mid)
    return hi


def anneal(D, tour, target, rng, end_temperature):
    """Geometric cooling from the solved T0 to end_temperature over exactly BUDGET evaluated moves."""
    tour = list(tour)
    length = tour_length(D, tour)
    t0 = solve_temperature(sample_worsening(D, tour, rng), target)
    temperature, cool = t0, (end_temperature / t0) ** (1 / BUDGET) if t0 > end_temperature else 1.0
    for _ in range(BUDGET):
        i, j = sorted(rng.sample(range(1, CITIES), 2))
        delta = two_opt_delta(D, tour, i, j)
        if delta < 0 or rng.random() < math.exp(-delta / temperature):
            tour[i:j + 1] = tour[i:j + 1][::-1]
            length += delta
        temperature *= cool
    assert sorted(tour) == list(range(CITIES)), "invalid tour"          # report only valid completed tours
    assert abs(length - tour_length(D, tour)) < 1e-6, "tracked length drifted"
    return length, t0


def run(target_acceptance):
    ratio = {"random": [], "nn": []}
    t0_scaled = {"random": [], "nn": []}
    naive_acceptance = []
    for k in range(INSTANCES):
        D = make_instance(SEED * 100 + k)
        rng = random.Random(SEED * 100 + k)
        nn = nearest_neighbour_tour(D)
        base = tour_length(D, nn)                                    # reference: the nearest-neighbour tour, no search
        edge = base / CITIES
        shuffled = list(range(CITIES))
        rng.shuffle(shuffled)
        for name, start in (("random", shuffled), ("nn", nn)):
            length, t0 = anneal(D, start, target_acceptance, rng, END_TEMPERATURE * edge)
            ratio[name].append(length / base)
            t0_scaled[name].append(t0 / edge)
        # the shortcut the chapter warns about: set T from the mean worsening delta so that exp(-mean/T) equals the target
        deltas = sample_worsening(D, nn, rng)
        naive_t = -(sum(deltas) / len(deltas)) / math.log(target_acceptance)
        naive_acceptance.append(acceptance(deltas, naive_t))
    return clean({
        "target_acceptance": target_acceptance,
        "ratio": {k: float(np.mean(v)) for k, v in ratio.items()},
        "per_instance": ratio,
        "t0_over_edge": {k: float(np.mean(v)) for k, v in t0_scaled.items()},
        "naive_acceptance": float(np.mean(naive_acceptance)),
        "cities": CITIES, "budget": BUDGET, "instances": INSTANCES,
    })

import pandas as pd
from kaggle_companion.activities.ch54 import explain

VALUES = [0.01, 0.1, 0.5, 0.9]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch54 import draw, NCOLS, HEIGHT

DEFAULT = 0.5
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()