Problems and Algorithms Lab

# Counting search spaces (Section 2.1.1)
import math

# Fitness plan: choose r = 3 of n = 5 exercises
n, r = 5, 3
print("combinations with repetition:   ", math.comb(n + r - 1, r))       # 35
print("combinations without repetition:", math.comb(n, r))               # 10, "5 choose 3"
print("permutations with repetition:   ", n ** r)                        # 125
print("permutations without repetition:", math.perm(n, r))               # 60

# Traveling salesman problem: brute force vs. dynamic programming (Bellman-Held-Karp)
for cities in (5, 10, 13, 20):
    print(f"{cities:2d} cities: {math.factorial(cities):.3e} tours (ATSP), "
          f"{math.factorial(cities - 1) // 2:.3e} distinct symmetric tours, Held-Karp ~ {2 ** cities * cities ** 2:.3e} steps")


def stirling2(n, k):
    """Stirling number of the second kind: ways to partition n objects into k non-empty groups (equation 2.1)."""
    return sum((-1) ** i * math.comb(k, i) * (k - i) ** n for i in range(k + 1)) // math.factorial(k)


# Smart cart clustering: n carts around k assembly points
print(f"S(50, 4) = {stirling2(50, 4):.2e}")                              # 5.3e28
print(f"S(100, 4) = {stirling2(100, 4):.2e}")                            # 6.7e58

# Unimodal, multimodal and deceptive landscapes; deterministic vs. stochastic search (Listing 2.1, Section 2.2)
import math
import random


def t(a):
    """|sin(pi a)| / (pi a), the factor in Listing 2.1, taken as 1 at a = 0 (the formula divides by zero there)."""
    return 1.0 if abs(a) < 1e-12 else abs(math.sin(math.pi * a)) / (math.pi * a)


unimodal = lambda x, y: x ** 2 + y ** 2
multimodal = lambda x, y: math.sin(x) * math.cos(y)
deceptive = lambda x, y: (1 - (t(x - 2) * t(y - 2)) ** 5) * (2 + (x - 7) ** 2 + 2 * (y - 7) ** 2)   # Listing 2.1


def compass_search(f, x, y, step=0.5, lo=-3, hi=3):
    """Deterministic local search: try the four compass moves, halve the step when none improves."""
    moves = 0
    while step > 1e-4:
        cand = [(x + dx, y + dy) for dx, dy in ((step, 0), (-step, 0), (0, step), (0, -step)) if lo <= x + dx <= hi and lo <= y + dy <= hi]
        best = min(cand, key=lambda p: f(*p))
        if f(*best) < f(x, y):
            (x, y), moves = best, moves + 1
        else:
            step /= 2
    return x, y


def random_walk_search(f, x, y, rng, sigma=0.4, iters=400, lo=-3, hi=3):
    """Stochastic local search: accept a random nearby point when it is better."""
    for _ in range(iters):
        nx, ny = min(hi, max(lo, x + rng.gauss(0, sigma))), min(hi, max(lo, y + rng.gauss(0, sigma)))
        if f(nx, ny) < f(x, y):
            x, y = nx, ny
    return x, y


for name, f in (("unimodal", unimodal), ("multimodal", multimodal), ("deceptive", deceptive)):
    a, b = compass_search(f, 2.5, 2.5), compass_search(f, 2.5, 2.5)
    print(f"{name:10s} deterministic, two runs: ({a[0]:.3f}, {a[1]:.3f}) and ({b[0]:.3f}, {b[1]:.3f}), f = {f(*a):.4f}")
    for seed in (1, 2):
        x, y = random_walk_search(f, 2.5, 2.5, random.Random(seed))
        print(f"{'':10s} stochastic, seed {seed}: ({x:.3f}, {y:.3f}), f = {f(x, y):.4f}")

# The n-queen constraint-satisfaction problem (Section 2.1.2, Listing 2.2, Figure 2.5)
import random


def attacks(q, c):
    """Columns whose queens attack the queen in column c (same row or same diagonal)."""
    return [j for j in range(len(q)) if j != c and (q[j] == q[c] or abs(q[j] - q[c]) == abs(j - c))]


def conflicts(q):
    """q[c] is the row of the queen in column c; count attacking pairs."""
    return sum(len(attacks(q, c)) for c in range(len(q))) // 2


def book_sliding(n, max_steps=500):
    """Listing 2.2: place column by column; when a column is impossible, drop the first column and slide the board."""
    board, col, seen = [0] + [None] * (n - 1), 1, set()
    for step in range(1, max_steps):
        if col == n:
            return board, step
        if (tuple(board), col) in seen:
            return None, step                          # the same state again: the procedure cycles forever
        seen.add((tuple(board), col))
        placed = [(c, board[c]) for c in range(col)]
        for row in range(n):
            if all(row != r and abs(row - r) != col - c for c, r in placed):
                board[col] = row
                col += 1
                break
        else:
            board = board[1:] + [None]
            col -= 1
    return None, max_steps


def min_conflicts(n, rng, max_moves=1000):
    """Heuristic repair (Figure 2.5): move a conflicted queen to the row with the fewest conflicts."""
    q = [rng.randrange(n) for _ in range(n)]
    for move in range(max_moves):
        bad = [c for c in range(n) if attacks(q, c)]
        if not bad:
            return q, move
        c = rng.choice(bad)
        q[c] = min(range(n), key=lambda r: (len(attacks(q[:c] + [r] + q[c + 1:], c)), rng.random()))
    return None, max_moves


def count_solutions(n, col=0, rows=(), d1=(), d2=()):
    """Backtracking search that counts every solution."""
    if col == n:
        return 1
    return sum(count_solutions(n, col + 1, rows + (r,), d1 + (r - col,), d2 + (r + col,))
               for r in range(n) if r not in rows and r - col not in d1 and r + col not in d2)


for n in range(4, 9):
    sol, steps = book_sliding(n)
    q, k = min_conflicts(n, random.Random(n))
    print(f"n = {n}: Listing 2.2 {'reaches ' + str(sol) + ' in ' + str(steps) + ' steps' if sol else 'cycles and never finishes'}; "
          f"min-conflicts needs {k} moves; {count_solutions(n)} solutions exist")

# Constrained optimization: Lagrange multipliers and linear programming (Sections 2.1.3 and 2.1.4)
from itertools import combinations

# 1. Ticket sales with an equality constraint x + y = 10,000 (Listing 2.5)
f = lambda x, y: 155 * x + (0.001 * x ** 1.5 + 70) * y - 1_000_000
g = lambda x: f(x, 10_000 - x)                  # substitute the constraint: one variable is left
lo, hi = 0.0, 10_000.0
for _ in range(200):                            # ternary search on the profit along the constraint line
    m1, m2 = lo + (hi - lo) / 3, hi - (hi - lo) / 3
    lo, hi = (m1, hi) if g(m1) < g(m2) else (lo, m2)
x = (lo + hi) / 2
y = 10_000 - x
lam = 155 + 0.0015 * x ** 0.5 * y               # dL/dx = 0 gives lambda = df/dx
print(f"in person x = {x:.0f}, online y = {y:.0f}, profit = ${f(x, y):,.0f}, lambda = {lam:.2f} (= df/dy = {0.001 * x ** 1.5 + 70:.2f})")

# 2. Animal feed mix LP (Table 2.1, Listing 2.6). With limestone = 1 - corn - soy, the problem lives in the
#    (corn, soy) plane; each constraint is written as a*c + b*s <= r.
cons = {
    "calcium >= 0.8%": (0.379, 0.378, 0.38 - 0.008),   # calcium = 0.38 - 0.379c - 0.378s
    "calcium <= 1.2%": (-0.379, -0.378, 0.012 - 0.38),
    "protein >= 22%": (-0.09, -0.5, -0.22),
    "fiber <= 5%": (0.02, 0.08, 0.05),
    "corn >= 0": (-1, 0, 0),
    "soybean meal >= 0": (0, -1, 0),
    "limestone >= 0": (1, 1, 1),
}
cost = lambda c, s: 30.5 * c + 10 * (1 - c - s) + 90 * s    # cents per kg
feasible = lambda c, s: all(a * c + b * s <= r + 1e-9 for a, b, r in cons.values())
vertices = []
for (a1, b1, r1), (a2, b2, r2) in combinations(cons.values(), 2):
    det = a1 * b2 - a2 * b1
    if abs(det) > 1e-12:
        c, s = (r1 * b2 - r2 * b1) / det, (a1 * r2 - a2 * r1) / det
        if feasible(c, s):
            vertices.append((c, s))
c, s = min(vertices, key=lambda v: cost(*v))    # a linear objective is optimal at a vertex of the feasible region
print(f"{len(vertices)} feasible vertices; best: corn {c:.2%}, limestone {1 - c - s:.2%}, soybean meal {s:.2%}, "
      f"cost ${cost(c, s) / 100:.4f} per kg")

# Heuristics: bin packing and the cargo bike search strategies (Section 2.3)

# 1. One-dimensional bin packing (Figure 2.17)
def first_fit(items, cap=7):
    bins = []
    for it in items:
        for b in bins:
            if sum(b) + it <= cap:
                b.append(it)
                break
        else:
            bins.append([it])
    return bins


items = [2, 4, 3, 6, 5, 1]
print("lower bound:", -(-sum(items) // 7), "bins")
print("first fit:           ", first_fit(items))
print("first fit decreasing:", first_fit(sorted(items, reverse=True)))

# 2. Loading a 4 x 8 container (Listing 2.7, written without NumPy)
W, H = 4, 8


def load(objects):
    grid = [[0] * W for _ in range(H)]

    def fit(w, h, k, rotate=True):
        for i in range(H - h + 1):              # move down first, then right
            for j in range(W - w + 1):
                if all(grid[i + a][j + b] == 0 for a in range(h) for b in range(w)):
                    for a in range(h):
                        for b in range(w):
                            grid[i + a][j + b] = k
                    return True
        return rotate and fit(h, w, k, False)   # last resort: rotate by 90 degrees

    left_out = [o for k, o in enumerate(objects, 1) if not fit(*o, k)]
    return grid, left_out


objects = [[3, 1], [3, 3], [5, 1], [4, 2], [3, 2]]   # [width, height]
for name, order in (("largest first", sorted(objects, key=lambda o: o[0] * o[1], reverse=True)),
                    ("smallest width first", sorted(objects, key=lambda o: (o[0], o[0] * o[1])))):
    grid, left_out = load(order)
    print(f"{name}: left out {left_out}")
    for row in reversed(grid):
        print("   ", " ".join(str(v) if v else "." for v in row))

# 3. Cargo bike: construction, modification and recombination (Tables 2.2 and 2.3, Figure 2.21)
weight = [14.6, 20, 8.5, 10, 13, 9.6, 4.9, 16.5, 8.77, 7.8]
profit = [14.54, 15.26, 5.8, 12.12, 8.2, 7.4, 10.3, 13.5, 6.6, 20.9]
total = lambda x: (round(sum(w for w, b in zip(weight, x) if b), 2), round(sum(p for p, b in zip(profit, x) if b), 2))

x, load_kg = [0] * 10, 0.0
for i in sorted(range(10), key=lambda i: profit[i] / weight[i], reverse=True):    # greedy by efficiency ($/kg)
    if load_kg + weight[i] <= 100:
        x[i] = 1
        load_kg += weight[i]
print("construction: ", tuple(x), "kg, $ =", total(x))
y = x[:]
y[8], y[4] = 0, 1                                                              # take item 9 out, put item 5 in
print("modification: ", tuple(y), "kg, $ =", total(y))
S1, S2 = (1, 1, 1, 1, 1, 0, 0, 1, 0, 1), (0, 1, 0, 1, 1, 0, 1, 1, 1, 1)
S3 = S1[:2] + S2[2:]                                                           # first two genes of S1, last eight of S2
print("recombination:", S1, total(S1), "+", S2, total(S2), "->", S3, total(S3))