NP-Completeness & Approximation Algorithms

# ── 3SAT VERIFIER ── O(clauses * 3) ──────────────────────────────────
def verify_3sat(clauses, assignment):
    """
    Verify a 3SAT certificate in polynomial time.
    clauses: list of 3-tuples like [(1, -2, 3), (-1, 4, -5)]
             positive = variable, negative = negated variable.
    assignment: dict {var: True/False}.
    O(clauses * 3) = O(n): clearly polynomial.
    """
    for clause in clauses:
        # Each clause is satisfied if at least one literal is true
        if not any(
            (assignment.get(abs(lit), False) if lit > 0
             else not assignment.get(abs(lit), False))
            for lit in clause
        ):
            return False   # clause not satisfied
    return True

# Example: (x1 OR x3 OR NOT x2) AND (NOT x1 OR x4 OR NOT x5)
clauses    = [(1, 3, -2), (-1, 4, -5)]
assignment = {1: True, 2: False, 3: True, 4: True, 5: False}
print("3SAT verified:", verify_3sat(clauses, assignment))   # True


# ── SUBSET SUM (weak NP) ── pseudo-polynomial O(n*T) DP ──────────────
def subset_sum_dp(nums, target):
    """
    Pseudo-polynomial DP: O(n*T) time, O(T) space.
    Works because T is treated as a number, not its bit-length.
    If T = 2^n, this is exponential in the input size (weak NP-hardness).
    """
    dp = {0}    # set of achievable sums
    for x in nums:
        dp = dp | {s + x for s in dp if s + x <= target}
    return target in dp

print("Subset sum:", subset_sum_dp([3,1,4,1,5,9,2,6], 11))  # True


# ── VERTEX COVER VERIFIER ── O(E) ────────────────────────────────────
def verify_vertex_cover(graph_edges, cover):
    """Verify vertex cover: every edge has at least one endpoint in cover."""
    cover_set = set(cover)
    return all(u in cover_set or v in cover_set for u, v in graph_edges)

# ── REDUCTION: Vertex Cover → Set Cover ──────────────────────────────
def vc_to_set_cover(n_vertices, edges):
    """
    Karp reduction: Vertex Cover instance → Set Cover instance.
    Universe S = E (set of edges, indexed 0..m-1).
    Set_v = {edges incident to vertex v}.
    k-vertex cover ↔ k-set cover.
    Polynomial time: O(V + E).
    """
    universe = list(range(len(edges)))    # edges are elements
    sets = {}
    for v in range(n_vertices):
        sets[v] = [i for i, (u, w) in enumerate(edges) if u==v or w==v]
    return universe, sets

edges = [(0,1),(0,2),(1,2),(1,3),(2,4),(3,4)]
universe, sets = vc_to_set_cover(5, edges)
print(f"Universe (edges): {universe}")
print(f"Sets (per vertex): {sets}")

import heapq, math
from collections import defaultdict

# ── 2-APPROX VERTEX COVER ─────────────────────────────────────────────
def approx_vertex_cover_2(edges):
    """
    2-approximation for Minimum Vertex Cover.
    Pick any edge, add both endpoints, remove incident edges.
    C <= 2 * C_opt. Proven tight (matching lower bound exists).
    O(E) time.
    """
    cover     = set()
    remaining = list(edges)
    covered   = set()

    for u, v in remaining:
        if (u, v) in covered or (v, u) in covered:
            continue                    # edge already covered
        cover.add(u)
        cover.add(v)
        # Mark all edges incident to u or v as covered
        for a, b in remaining:
            if a == u or b == u or a == v or b == v:
                covered.add((a, b))
    return cover


# ── GREEDY SET COVER (ln n + 1 approximation) ─────────────────────────
def greedy_set_cover(universe, sets):
    """
    Greedy set cover: always pick the set covering the most uncovered elements.
    Approximation ratio: ln(|universe|) + 1.
    O(|universe| * |sets|) per iteration → O(n^2 * m) total.
    Lower bound: no poly-time alg achieves (1-eps)*ln(n) unless NP ⊆ DTIME(n^(log log n)).
    """
    uncovered = set(universe)
    cover     = []

    while uncovered:
        # Pick set with maximum intersection with uncovered elements
        best = max(sets, key=lambda s: len(uncovered & sets[s]))
        cover.append(best)
        uncovered -= sets[best]

    return cover


# ── METRIC TSP: MST-BASED 2-APPROXIMATION ─────────────────────────────
def metric_tsp_2approx(n, dist):
    """
    2-approximation for Metric TSP via MST + DFS preorder.
    dist[i][j] must satisfy triangle inequality.
    Returns (tour_cost, tour_order).
    O(n^2) Prim + O(n) DFS = O(n^2) total.
    """
    # Step 1: Prim's MST (O(n^2) with adjacency matrix)
    in_mst = [False]*n
    key    = [float('inf')]*n
    parent = [-1]*n
    key[0] = 0
    adj    = defaultdict(list)  # MST adjacency list

    for _ in range(n):
        # Find min-key vertex not in MST
        u = min((i for i in range(n) if not in_mst[i]), key=lambda i: key[i])
        in_mst[u] = True
        if parent[u] != -1:
            adj[parent[u]].append(u)
            adj[u].append(parent[u])
        for v in range(n):
            if not in_mst[v] and dist[u][v] < key[v]:
                key[v] = dist[u][v]; parent[v] = u

    # Step 2: DFS preorder traversal (shortcutting = visiting unvisited)
    visited = [False]*n
    tour    = []

    def dfs(u):
        visited[u] = True
        tour.append(u)
        for v in sorted(adj[u]):   # sorted for determinism
            if not visited[v]:
                dfs(v)

    dfs(0)

    # Step 3: Compute tour cost (triangle inequality shortcutting is implicit)
    tour_cost = sum(dist[tour[i]][tour[i+1]] for i in range(len(tour)-1))
    tour_cost += dist[tour[-1]][tour[0]]   # return to start

    return tour_cost, tour

def partition_ptas(items, epsilon):
    """
    PTAS for load-balancing partition: minimise max(w(A), w(B)).
    Phase 1: solve exactly for largest m = ceil(1/epsilon) - 1 items.
    Phase 2: greedy assign remaining items.
    Approximation ratio: 1 + epsilon.
    Time: O(2^(1/epsilon) + n) — polynomial in n, exponential in 1/epsilon.
    """
    items_sorted = sorted(items, reverse=True)
    m = max(1, math.ceil(1/epsilon) - 1)

    # Phase 1: brute-force optimal partition of first m items
    best_diff = float('inf')
    best_mask = 0
    for mask in range(1 << min(m, len(items_sorted))):
        wA = sum(items_sorted[i] for i in range(min(m, len(items_sorted)))
                 if mask & (1 << i))
        wB = sum(items_sorted[i] for i in range(min(m, len(items_sorted)))
                 if not mask & (1 << i))
        if abs(wA - wB) < best_diff:
            best_diff = abs(wA - wB)
            best_mask = mask

    A = [items_sorted[i] for i in range(min(m, len(items_sorted)))
         if best_mask & (1 << i)]
    B = [items_sorted[i] for i in range(min(m, len(items_sorted)))
         if not best_mask & (1 << i)]

    # Phase 2: greedy assignment of remaining items
    for item in items_sorted[m:]:
        if sum(A) <= sum(B):
            A.append(item)
        else:
            B.append(item)

    return max(sum(A), sum(B)), A, B


def subset_sum_fptas(nums, target, epsilon):
    """
    FPTAS for Subset Sum approximation: find subset with sum in
    [(1-epsilon)*target, target] if it exists.
    Time: O(n^2 / epsilon) — polynomial in both n and 1/epsilon.
    Strategy: round down all numbers by factor K=floor(epsilon*max/n),
              run exact DP on rounded numbers, unscale answer.
    """
    if not nums: return 0, []
    max_val = max(nums)
    K = max(1, int(epsilon * max_val / len(nums)))
    rounded = [x // K for x in nums]
    target_r = target // K

    # Exact DP on rounded numbers — O(n * target/K) = O(n^2/epsilon)
    dp = {0: []}
    for i, x in enumerate(rounded):
        new_dp = {}
        for s, path in dp.items():
            ns = s + x
            if ns <= target_r and ns not in new_dp:
                new_dp[ns] = path + [i]
        dp.update(new_dp)

    # Find best achievable sum
    best = max((s for s in dp if s*K <= target), default=0)
    chosen = [nums[i] for i in dp[best]]
    return sum(chosen), chosen