Blind Search Lab

# Blind search on the 8-puzzle (Chapter 3, Section 3.3)
# A state is a tuple of 9 tiles read row by row; 0 is the blank. Moves slide the blank left, up, right or down.
from collections import deque
import heapq

MOVES = [("Left", -1), ("Up", -3), ("Right", 1), ("Down", 3)]      # the operator order used in the book


def expand(s):
    i = s.index(0)
    for name, d in MOVES:
        j = i + d
        if name == "Left" and i % 3 == 0 or name == "Right" and i % 3 == 2 or not 0 <= j < 9:
            continue
        t = list(s)
        t[i], t[j] = t[j], t[i]
        yield name, tuple(t)


def solvable(s):
    tiles = [t for t in s if t]
    return sum(a > b for k, a in enumerate(tiles) for b in tiles[k + 1:]) % 2 == 0


def path_to(parent, s):
    p = []
    while s is not None:
        p.append(s)
        s = parent[s]
    return p[::-1]


def graph_search(start, goal, lifo=False):
    """BFS (FIFO queue) or DFS (LIFO stack), with the goal test on generation as in algorithms 3.1 and 3.2."""
    frontier, parent, explored = deque([start]), {start: None}, set()
    in_frontier = {start}
    while frontier:
        s = frontier.pop() if lifo else frontier.popleft()
        in_frontier.discard(s)
        explored.add(s)
        for _, c in expand(s):
            if c not in explored and c not in in_frontier:
                parent[c] = s
                if c == goal:
                    return path_to(parent, c), len(explored)
                frontier.append(c)
                in_frontier.add(c)
    return None, len(explored)


def depth_limited(start, goal, limit):
    """Recursive DLS that never moves back to a state already on the current path. Returns (path, visits)."""
    visits = 0

    def rec(s, path):
        nonlocal visits
        visits += 1
        if s == goal:
            return path
        if len(path) - 1 == limit:
            return None
        for _, c in expand(s):
            if c not in path:
                r = rec(c, path + [c])
                if r:
                    return r
        return None

    return rec(start, [start]), visits


def iterative_deepening(start, goal, max_depth=30):
    total = 0
    for limit in range(max_depth + 1):
        p, v = depth_limited(start, goal, limit)
        total += v
        if p:
            return p, total, limit
    return None, total, None


def ucs(start, goal):
    """Uniform-cost search: every move costs 1, so it expands states in the same order as BFS."""
    pq, best, parent, popped, k = [(0, 0, start)], {start: 0}, {start: None}, set(), 1
    while pq:
        g, _, s = heapq.heappop(pq)
        if s in popped:
            continue
        popped.add(s)
        if s == goal:
            return path_to(parent, s), len(popped)
        for _, c in expand(s):
            if c not in popped and g + 1 < best.get(c, 1e9):
                best[c], parent[c] = g + 1, s
                heapq.heappush(pq, (g + 1, k, c)); k += 1
    return None, len(popped)


start = (1, 4, 2, 3, 7, 5, 6, 0, 8)       # Listing 3.9
goal = (0, 1, 2, 3, 4, 5, 6, 7, 8)
print("solvable:", solvable(start))
p, n = graph_search(start, goal)
print(f"BFS: {len(p) - 1} moves, {n} states explored")
p, n = graph_search(start, goal, lifo=True)
print(f"DFS: {len(p) - 1} moves, {n} states explored")
p, v = depth_limited(start, goal, 5)
print(f"DLS (limit 5): {len(p) - 1} moves, {v} states visited")
p, v, L = iterative_deepening(start, goal)
print(f"IDS: {len(p) - 1} moves, {v} states visited, depth {L}")
p, n = ucs(start, goal)
print(f"UCS: {len(p) - 1} moves, {n} states explored")
for s in p:
    print(s[:3], s[3:6], s[6:], sep="\n", end="\n\n")

# BFS and DFS path planning for a mobile robot on a grid with obstacles (Figures 3.13 to 3.25)
from collections import deque
import matplotlib.pyplot as plt

GRID = ["S.#.....",
        ".##.##..",
        "....#...",
        ".##...#.",
        "...#.#..",
        ".#...#.G"]                     # S start, G destination, # obstacle
R, C = len(GRID), len(GRID[0])
start = next((r, c) for r in range(R) for c in range(C) if GRID[r][c] == "S")
goal = next((r, c) for r in range(R) for c in range(C) if GRID[r][c] == "G")
# compass moves, tried in this order: S, SE, E, NE, N, NW, W, SW
MOVES = [(1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1)]


def neighbors(cell):
    r, c = cell
    for dr, dc in MOVES:
        rr, cc = r + dr, c + dc
        if 0 <= rr < R and 0 <= cc < C and GRID[rr][cc] != "#":
            yield rr, cc


def search(lifo):
    frontier, parent, explored = deque([start]), {start: None}, []
    while frontier:
        cell = frontier.pop() if lifo else frontier.popleft()   # LIFO stack for DFS, FIFO queue for BFS
        explored.append(cell)
        for n in neighbors(cell):
            if n not in parent:
                parent[n] = cell
                if n == goal:
                    path = [n]
                    while parent[path[-1]]:
                        path.append(parent[path[-1]])
                    return path[::-1], explored
                frontier.append(n)
    return None, explored


fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
for ax, (name, lifo) in zip(axes, [("BFS", False), ("DFS", True)]):
    path, explored = search(lifo)
    ax.imshow([[1 if ch == "#" else 0 for ch in row] for row in GRID], cmap="Greys")
    ax.scatter([c for r, c in explored], [r for r, c in explored], c="tab:orange", s=40, label="explored")
    ax.plot([c for r, c in path], [r for r, c in path], "b-o", label="path")
    ax.set_title(f"{name}: {len(path) - 1} moves, {len(explored)} cells explored")
    ax.legend(loc="lower left")
    print(name, "path:", path)
plt.show()

# Dijkstra, uniform-cost search and bidirectional Dijkstra on the example of Figures 3.27 to 3.32
import heapq
import math

EDGES = [("S", "A", 11), ("S", "C", 3), ("S", "D", 6), ("C", "A", 6), ("C", "B", 8), ("C", "D", 2), ("A", "B", 3), ("D", "B", 5)]
G = {}
for u, v, w in EDGES:                       # undirected graph
    G.setdefault(u, []).append((v, w))
    G.setdefault(v, []).append((u, w))


def dijkstra(G, source):
    """Original Dijkstra: shortest distances from the source to every node (algorithm 3.3)."""
    dist = {n: math.inf for n in G}
    pred = {n: None for n in G}
    dist[source] = 0
    unrelaxed = set(G)
    while unrelaxed:
        node = min(unrelaxed, key=dist.get)                 # the closest unrelaxed node
        unrelaxed.remove(node)
        for child, w in G[node]:
            if child in unrelaxed and dist[node] + w < dist[child]:
                dist[child], pred[child] = dist[node] + w, node
        print(f"relax {node}: " + ", ".join(f"{n}={dist[n]}({pred[n] or '-'})" for n in sorted(G)))
    return dist, pred


def route(pred, target):
    r = []
    while target:
        r.append(target)
        target = pred[target]
    return r[::-1]


dist, pred = dijkstra(G, "S")
print("S to A:", "-".join(route(pred, "A")), dist["A"])           # S-C-A, 9
print("S to B:", "-".join(route(pred, "B")), dist["B"])           # S-C-D-B, 10


def ucs(G, source, goal):
    """Uniform-cost search: fills the priority queue as it goes and stops when the goal is popped."""
    pq, best, pred, popped = [(0, source)], {source: 0}, {source: None}, []
    while pq:
        g, node = heapq.heappop(pq)
        if node in popped:
            continue
        popped.append(node)
        if node == goal:
            return route(pred, goal), g, popped
        for child, w in G[node]:
            if g + w < best.get(child, math.inf):
                best[child], pred[child] = g + w, node
                heapq.heappush(pq, (g + w, child))


r, c, popped = ucs(G, "S", "B")
print("UCS S to B:", "-".join(r), c, " popped in order:", popped)


def bidirectional_dijkstra(G, s, t):
    """Search forward from s and backward from t; stop when the two smallest keys add up to at least the best meeting cost."""
    dist = [{s: 0}, {t: 0}]
    pred = [{s: None}, {t: None}]
    pq = [[(0, s)], [(0, t)]]
    done = [set(), set()]
    best, meet, side = math.inf, None, 0
    while pq[0] and pq[1]:
        if pq[0][0][0] + pq[1][0][0] >= best:                # no shorter path can still be found
            break
        g, u = heapq.heappop(pq[side])
        if u not in done[side]:
            done[side].add(u)
            for v, w in G[u]:
                if g + w < dist[side].get(v, math.inf):
                    dist[side][v], pred[side][v] = g + w, u
                    heapq.heappush(pq[side], (g + w, v))
                if v in dist[1 - side] and g + w + dist[1 - side][v] < best:
                    best, meet = g + w + dist[1 - side][v], v
        side = 1 - side                                      # alternate directions
    forward = route(pred[0], meet)
    backward = route(pred[1], meet)[::-1][1:]
    return forward + backward, best


print("Bidirectional S to B:", bidirectional_dijkstra(G, "S", "B"))

# Blind search on a real road network: University of Toronto, St. George campus (Section 3.5)
# !pip install osmnx          # uncomment in Colab
from collections import deque
import heapq
import math
import osmnx

G = osmnx.graph_from_point((43.661667, -79.395), dist=300, simplify=True)   # around King's College Circle
origin = osmnx.distance.nearest_nodes(G, -79.392442, 43.664527)                # King Edward VII equestrian statue
destination = osmnx.distance.nearest_nodes(G, -79.397669, 43.659659)           # Bahen Centre
# keep the shortest of any parallel edges and ignore self-loops (Figures 3.38 and 3.39)
adj = {n: {} for n in G.nodes}
for u, v, data in G.edges(data=True):
    if u != v and data["length"] < adj[u].get(v, math.inf):
        adj[u][v] = data["length"]
radj = {n: {} for n in G.nodes}
for u in adj:
    for v, w in adj[u].items():
        radj[v][u] = w


def route_of(parent, n):
    r = []
    while n is not None:
        r.append(n)
        n = parent[n]
    return r[::-1]


def length(r):
    return sum(adj[a][b] for a, b in zip(r, r[1:]))


def traverse(lifo):
    """BFS (queue) or DFS (stack): find a path, ignoring edge lengths."""
    frontier, parent, explored = deque([origin]), {origin: None}, set()
    while frontier:
        n = frontier.pop() if lifo else frontier.popleft()
        explored.add(n)
        for c in adj[n]:
            if c not in parent:
                parent[c] = n
                if c == destination:
                    return route_of(parent, c), len(explored)
                frontier.append(c)


def ucs():
    """Uniform-cost search (Dijkstra that stops at the destination)."""
    pq, best, parent, done = [(0, origin)], {origin: 0}, {origin: None}, set()
    while pq:
        g, n = heapq.heappop(pq)
        if n in done:
            continue
        done.add(n)
        if n == destination:
            return route_of(parent, n), len(done)
        for c, w in adj[n].items():
            if g + w < best.get(c, math.inf):
                best[c], parent[c] = g + w, n
                heapq.heappush(pq, (g + w, c))


def bidirectional():
    """Bidirectional Dijkstra with the correct stopping rule: stop when top_f + top_b >= best meeting cost."""
    graphs, dist, parent = (adj, radj), [{origin: 0}, {destination: 0}], [{origin: None}, {destination: None}]
    pq, done, best, meet, side = [[(0, origin)], [(0, destination)]], [set(), set()], math.inf, None, 0
    while pq[0] and pq[1] and pq[0][0][0] + pq[1][0][0] < best:
        g, n = heapq.heappop(pq[side])
        if n not in done[side]:
            done[side].add(n)
            for c, w in graphs[side][n].items():
                if g + w < dist[side].get(c, math.inf):
                    dist[side][c], parent[side][c] = g + w, n
                    heapq.heappush(pq[side], (g + w, c))
                if c in dist[1 - side] and g + w + dist[1 - side][c] < best:
                    best, meet = g + w + dist[1 - side][c], c
        side = 1 - side
    r = route_of(parent[0], meet) + route_of(parent[1], meet)[::-1][1:]
    return r, len(done[0] | done[1])


for name, fn in [("BFS", lambda: traverse(False)), ("DFS", lambda: traverse(True)), ("UCS / Dijkstra", ucs), ("Bidirectional Dijkstra", bidirectional)]:
    r, n = fn()
    print(f"{name:23s} {length(r):8.1f} m, {n:4d} nodes explored")
osmnx.plot_graph_route(G, ucs()[0], route_color="r", node_size=0)

# Running time by complexity class (Table 3.3) and BFS time and memory (Table 3.4)
import math

n, ops = 20_000, 1e6                                  # problem size, operations per second
for name, f in [("O(1)", lambda n: 1), ("O(log n)", lambda n: math.log2(n)), ("O(n)", lambda n: n),
                ("O(n log n)", lambda n: n * math.log2(n)), ("O(n^2)", lambda n: n**2), ("O(n^3)", lambda n: n**3)]:
    s = f(n) / ops
    print(f"{name:11s} {s:.3g} s" + (f" = {s / 60:.2f} min" if s > 60 else "") + (f" = {s / 86400:.1f} days" if s > 86400 else ""))
print(f"O(2^n)      about 10^{n * math.log10(2) - math.log10(ops * 3.156e7):.0f} years")
print(f"O(n!)       about 10^{math.lgamma(n + 1) / math.log(10) - math.log10(ops * 3.156e7):.0f} years")

# BFS with branching factor b = 5, 1 KB per node, 1,000 nodes per second: N = (b^(d+1) - 1) / (b - 1)
b = 5
for d in range(2, 15, 2):
    N = (b ** (d + 1) - 1) // (b - 1)
    t = N / 1000
    when = f"{t:.3g} s" if t < 60 else f"{t / 60:.3g} min" if t < 3600 else f"{t / 3600:.3g} h" if t < 86400 else f"{t / 86400:.3g} days" if t < 3.156e7 else f"{t / 3.156e7:.3g} years"
    print(f"d = {d:2d}: N = {N:>15,d} nodes, time {when:>12s}, memory {N / 1e6:.3g} GB")