""" geo-optimize — neutral, product-agnostic optimization technique library. ARCHITECTURE (non-negotiable boundary) -------------------------------------- This is a *shared technique library*, NOT a shared runtime. It contains: - a placement / coverage / dispatch problem formulation, - a uniform solver interface (classical exact, classical heuristic, quantum-inspired, and staged real-QPU backends), - a benchmark harness that scores every backend with the SAME true objective. It contains NO product data. PHOENIX (Sicily wildfire), OLDGABE (US wildfire), and INGV (Etna science) each instantiate this library separately, supplying their own objectives, candidate sites, constraints, risk layers, and data planes. The instances never see each other's feeds. There is deliberately no function here that fuses two products' confidences — that would be a cross-contamination vector by design. Build the library once; customize per mission. HONEST STANCE ON QUANTUM ------------------------ At realistic problem sizes (dozens of candidate sites, hundreds of demand cells) maximum-coverage placement is solved to PROVABLE OPTIMALITY by CP-SAT in milliseconds. Quantum and quantum-inspired solvers are wired in as *swappable benchmark backends only*. We never report "quantum advantage" unless a backend beats the strongest classical baseline on the SAME instance by the SAME true objective. Documenting "no advantage at this scale" is itself the deliverable. COMPUTE ------- Everything here is CPU-light (tiny MILP + small annealing) and laptop-safe. The real-QPU path (D-Wave Advantage2) and any heavy sweep are staged for RunPod; see runpod/. Nothing in this module touches the DGX; if any benchmark is ever moved to DGX it must factor in load_guard per the standing rule. """ from __future__ import annotations import time import numpy as np from dataclasses import dataclass, field from typing import Callable # ========================================================================== # 1. Problem instance (the shared formulation) # ========================================================================== @dataclass class PlacementInstance: """ Maximum-coverage facility placement. candidates : C candidate sensor sites demand : D demand cells, each with a risk weight w_d >= 0 cover : (C, D) boolean — cover[c, d] is True if site c covers demand d (within detection range AND terrain line-of-sight) cost : (C,) install+maintenance cost proxy per site budget_p : max number of sites to place (cardinality budget) cost_cap : optional max total cost (None = ignore) allowed : (C,) boolean — site is eligible (e.g. LoRa-reachable to gateway, road-accessible). Ineligible sites are forced off. name : instance label (product/mission specific) meta : free-form dict (lat/lon, elevations, units, provenance) """ candidates: np.ndarray # (C, 2) local ENU metres, or arbitrary ids demand: np.ndarray # (D, 2) local ENU metres weights: np.ndarray # (D,) cover: np.ndarray # (C, D) bool cost: np.ndarray # (C,) budget_p: int cost_cap: float | None = None allowed: np.ndarray | None = None name: str = "instance" meta: dict = field(default_factory=dict) @property def C(self) -> int: return self.cover.shape[0] @property def D(self) -> int: return self.cover.shape[1] def eligible(self) -> np.ndarray: return np.ones(self.C, bool) if self.allowed is None else self.allowed.astype(bool) def total_risk(self) -> float: return float(self.weights.sum()) # ========================================================================== # 2. True objective (every backend is scored by THIS, regardless of search) # ========================================================================== def coverage_score(inst: PlacementInstance, selected) -> float: """Risk-weighted coverage of a chosen site set. The single source of truth.""" sel = np.asarray(list(selected), dtype=int) if sel.size == 0: return 0.0 covered = inst.cover[sel].any(axis=0) # (D,) bool return float(inst.weights[covered].sum()) def repair_to_feasible(inst: PlacementInstance, selected) -> list: """Project a raw (possibly infeasible) annealer decode onto the nearest feasible site set: drop ineligible sites, then greedily keep the budget_p highest-marginal sites within the cost cap. Standard post-processing for QUBO/QPU outputs.""" elig = inst.eligible() cand = [int(c) for c in selected if elig[c]] covered = np.zeros(inst.D, bool) chosen = [] while cand and len(chosen) < inst.budget_p: best_c, best_gain = None, -1.0 for c in cand: if inst.cost_cap is not None and \ inst.cost[chosen + [c]].sum() > inst.cost_cap + 1e-9: continue g = inst.weights[inst.cover[c] & ~covered].sum() if g > best_gain: best_gain, best_c = g, c if best_c is None: break chosen.append(best_c) covered |= inst.cover[best_c] cand.remove(best_c) return chosen def feasible(inst: PlacementInstance, selected) -> bool: sel = np.asarray(list(selected), dtype=int) if sel.size > inst.budget_p: return False if inst.cost_cap is not None and inst.cost[sel].sum() > inst.cost_cap + 1e-9: return False elig = inst.eligible() if sel.size and not elig[sel].all(): return False return True # ========================================================================== # 3. Solver result + uniform interface # ========================================================================== @dataclass class SolverResult: backend: str selected: list # chosen site indices score: float # true coverage_score latency_ms: float proven_optimal: bool = False # only exact solvers can assert this note: str = "" def _timed(fn, *a, **k): t0 = time.perf_counter() out = fn(*a, **k) return out, (time.perf_counter() - t0) * 1e3 # ---- 3a. Classical exact: OR-Tools CP-SAT (the deliverable) -------------- def solve_cpsat(inst: PlacementInstance, time_limit_s: float = 10.0) -> SolverResult: from ortools.sat.python import cp_model def _run(): m = cp_model.CpModel() elig = inst.eligible() x = [m.NewBoolVar(f"x{c}") for c in range(inst.C)] for c in range(inst.C): if not elig[c]: m.Add(x[c] == 0) y = [m.NewBoolVar(f"y{d}") for d in range(inst.D)] # demand d can be marked covered only if some covering site is chosen for d in range(inst.D): covering = [x[c] for c in range(inst.C) if inst.cover[c, d]] if covering: m.Add(y[d] <= sum(covering)) else: m.Add(y[d] == 0) m.Add(sum(x) <= inst.budget_p) if inst.cost_cap is not None: scale = 1000 m.Add(sum(int(round(inst.cost[c] * scale)) * x[c] for c in range(inst.C)) <= int(round(inst.cost_cap * scale))) wscale = 1000 m.Maximize(sum(int(round(inst.weights[d] * wscale)) * y[d] for d in range(inst.D))) solver = cp_model.CpSolver() solver.parameters.max_time_in_seconds = time_limit_s solver.parameters.num_search_workers = 8 status = solver.Solve(m) sel = [c for c in range(inst.C) if solver.Value(x[c]) == 1] optimal = status == cp_model.OPTIMAL return sel, optimal (sel, optimal), ms = _timed(_run) return SolverResult("cpsat_exact", sel, coverage_score(inst, sel), ms, proven_optimal=optimal, note="OR-Tools CP-SAT; provably optimal if proven_optimal=True") # ---- 3b. Classical heuristic: greedy (1 - 1/e) guarantee, deployable ------ def solve_greedy(inst: PlacementInstance) -> SolverResult: def _run(): elig = inst.eligible() covered = np.zeros(inst.D, bool) selected = [] cand = [c for c in range(inst.C) if elig[c]] while len(selected) < inst.budget_p and cand: best_c, best_gain = None, -1.0 for c in cand: if inst.cost_cap is not None and \ inst.cost[selected + [c]].sum() > inst.cost_cap + 1e-9: continue gain = inst.weights[inst.cover[c] & ~covered].sum() if gain > best_gain: best_gain, best_c = gain, c if best_c is None or best_gain <= 0: break selected.append(best_c) covered |= inst.cover[best_c] cand.remove(best_c) return selected sel, ms = _timed(_run) return SolverResult("greedy", sel, coverage_score(inst, sel), ms, note="classical greedy max-coverage; (1-1/e) approx; always-available fallback") # ---- 3c. Quantum-inspired: simulated annealing on the TRUE objective ----- def solve_sa_direct(inst: PlacementInstance, sweeps: int = 4000, restarts: int = 8, seed: int = 0) -> SolverResult: """Metaheuristic search over site subsets against the exact objective.""" def _run(): rng = np.random.default_rng(seed) elig = inst.eligible() elig_idx = np.where(elig)[0] best_sel, best_score = [], 0.0 for r in range(restarts): # start from greedy-ish random subset of size budget_p cur = set(rng.choice(elig_idx, size=min(inst.budget_p, elig_idx.size), replace=False).tolist()) cur_score = coverage_score(inst, cur) if feasible(inst, cur) else -1.0 T0, T1 = 1.0 * (inst.total_risk() / max(1, inst.D)), 1e-3 for s in range(sweeps): T = T0 * (T1 / T0) ** (s / max(1, sweeps - 1)) # propose: add / remove / swap move = rng.integers(0, 3) cand = set(cur) if move == 0 and len(cand) < inst.budget_p: cand.add(int(rng.choice(elig_idx))) elif move == 1 and cand: cand.discard(int(rng.choice(list(cand)))) elif cand: cand.discard(int(rng.choice(list(cand)))) cand.add(int(rng.choice(elig_idx))) if not feasible(inst, cand): continue cs = coverage_score(inst, cand) if cs >= cur_score or rng.random() < np.exp((cs - cur_score) / max(T, 1e-9)): cur, cur_score = cand, cs if cs > best_score: best_sel, best_score = sorted(cand), cs return best_sel sel, ms = _timed(_run) return SolverResult("sa_direct_qinspired", sel, coverage_score(inst, sel), ms, note="simulated annealing (quantum-inspired) on exact objective") # ---- 3d. Native-QUBO path: penalty QUBO + D-Wave Ocean SA (annealer proxy)- def build_coverage_qubo(inst: PlacementInstance, P_cov: float = None, P_bud: float = None): """ Max-coverage -> QUBO (minimize). Variables: x_c, y_d, plus slack bits to make the coverage inequality y_d <= sum_{c covers d} x_c exact. minimize -sum_d w_d y_d + P_cov * sum_d ( sum_{c in S_d} x_c - y_d - slack_d )^2 + P_bud * ( sum_c x_c - p )^2 (cardinality = p) Returns a dimod BinaryQuadraticModel and a decode() that maps a sample -> site set. This is the formulation a real D-Wave Advantage2 QPU consumes (see runpod/). """ import dimod w = inst.weights if P_cov is None: P_cov = 2.0 * float(w.max()) + 1.0 if P_bud is None: P_bud = 2.0 * float(w.sum()) + 1.0 bqm = dimod.BinaryQuadraticModel("BINARY") xv = [f"x{c}" for c in range(inst.C)] yv = [f"y{d}" for d in range(inst.D)] for v in xv + yv: bqm.add_variable(v, 0.0) # reward coverage for d in range(inst.D): bqm.add_variable(yv[d], -float(w[d])) # coverage inequality with slack bits per demand elig = inst.eligible() for d in range(inst.D): Sd = [c for c in range(inst.C) if inst.cover[c, d] and elig[c]] if not Sd: bqm.fix_variable(yv[d], 0) continue nslack = max(1, int(np.ceil(np.log2(len(Sd) + 1)))) slacks = [] for b in range(nslack): sv = f"s_{d}_{b}" bqm.add_variable(sv, 0.0) slacks.append((sv, 2 ** b)) # terms in g = (sum_x - y - sum slack) terms = [(xv[c], 1.0) for c in Sd] + [(yv[d], -1.0)] + \ [(sv, -coef) for sv, coef in slacks] # P_cov * g^2 expanded for i, (vi, ci) in enumerate(terms): bqm.add_variable(vi, P_cov * ci * ci) # ci^2 * vi (binary) for vj, cj in terms[i + 1:]: bqm.add_interaction(vi, vj, 2.0 * P_cov * ci * cj) # cardinality penalty P_bud (sum_x - p)^2 p = inst.budget_p for c in range(inst.C): bqm.add_variable(xv[c], P_bud * (1.0 - 2.0 * p)) for ci in range(inst.C): for cj in range(ci + 1, inst.C): bqm.add_interaction(xv[ci], xv[cj], 2.0 * P_bud) # force ineligible sites off for c in range(inst.C): if not elig[c]: bqm.fix_variable(xv[c], 0) def decode(sample): return sorted(c for c in range(inst.C) if xv[c] in sample and sample[xv[c]] == 1) return bqm, decode def solve_ocean_sa(inst: PlacementInstance, num_reads: int = 200, seed: int = 0) -> SolverResult: """D-Wave Ocean SimulatedAnnealingSampler on the native penalty QUBO. This is the canonical *classical proxy for the quantum annealer* — same QUBO the QPU would receive. Decoded x is scored by the TRUE objective.""" from dwave.samplers import SimulatedAnnealingSampler def _run(): bqm, decode = build_coverage_qubo(inst) ss = SimulatedAnnealingSampler().sample(bqm, num_reads=num_reads, seed=seed) # pick the decoded sample with best TRUE objective among feasible reads best_sel, best = [], -1.0 for rec in ss.data(["sample"]): sel = repair_to_feasible(inst, decode(rec.sample)) # decode + classical repair sc = coverage_score(inst, sel) if sc > best: best, best_sel = sc, sel return best_sel sel, ms = _timed(_run) return SolverResult("ocean_sa_qubo", sel, coverage_score(inst, sel), ms, note="D-Wave Ocean SA on native penalty-QUBO (annealer proxy)") # ---- 3e. Quantum-inspired: ballistic Simulated Bifurcation on the QUBO ---- def solve_sim_bifurcation(inst: PlacementInstance, steps: int = 600, restarts: int = 6, dt: float = 0.5, seed: int = 0, max_vars: int = 1200) -> SolverResult: """Ballistic Simulated Bifurcation (Goto et al.) — a distinct quantum-inspired solver. Runs on the same penalty-QUBO -> Ising. Gated by variable count to stay laptop-light; larger instances are skipped with a logged note.""" import dimod def _run(): bqm, decode = build_coverage_qubo(inst) if bqm.num_variables > max_vars: return None, f"skipped: {bqm.num_variables} vars > max_vars={max_vars} (stage on RunPod)" variables = list(bqm.variables) idx = {v: i for i, v in enumerate(variables)} n = len(variables) h = np.zeros(n) J = np.zeros((n, n)) lin, quad, off = bqm.to_ising() # dimod -> Ising: (h dict, J dict, offset) for v, b in lin.items(): h[idx[v]] = b for (u, v), c in quad.items(): J[idx[u], idx[v]] += c J[idx[v], idx[u]] += c rng = np.random.default_rng(seed) a0, c0 = 1.0, 0.5 / (np.sqrt(n) * np.std(J[J != 0]) + 1e-9) best_sel, best = [], -1.0 for r in range(restarts): x = 0.1 * (rng.random(n) - 0.5) y = 0.1 * (rng.random(n) - 0.5) for s in range(steps): a = a0 * s / steps # descend dimod Ising energy E(s)=h.s + s.J.s : force = -dE/dx y += dt * (-(a0 - a) * x - c0 * (J @ x + h)) x += dt * a0 * y over = np.abs(x) > 1.0 x[over] = np.sign(x[over]) y[over] = 0.0 spins = np.sign(x); spins[spins == 0] = 1 sample = {v: int((1 - spins[idx[v]]) / 2) for v in variables} # spin->binary sel = repair_to_feasible(inst, decode(sample)) # decode + classical repair sc = coverage_score(inst, sel) if sc > best: best, best_sel = sc, sel return best_sel, "" (out, note), ms = _timed(_run) if out is None: return SolverResult("sim_bifurcation_qinspired", [], 0.0, ms, note=note) return SolverResult("sim_bifurcation_qinspired", out, coverage_score(inst, out), ms, note="ballistic Simulated Bifurcation (quantum-inspired) on QUBO") # ========================================================================== # 4. Benchmark harness (shared; identical scoring for every backend) # ========================================================================== def benchmark(inst: PlacementInstance, backends: list[Callable] = None) -> dict: if backends is None: backends = [solve_cpsat, solve_greedy, solve_sa_direct, solve_ocean_sa, solve_sim_bifurcation] results: list[SolverResult] = [] for fn in backends: try: results.append(fn(inst)) except Exception as e: # a failed backend must not sink the run results.append(SolverResult(fn.__name__, [], 0.0, 0.0, note=f"ERROR: {e}")) # reference optimum = best proven-optimal score, else best observed exact = [r for r in results if r.proven_optimal] ref = max((r.score for r in exact), default=max(r.score for r in results)) quantum_like = {"sa_direct_qinspired", "ocean_sa_qubo", "sim_bifurcation_qinspired"} classical = {"cpsat_exact", "greedy"} best_classical = max((r.score for r in results if r.backend in classical), default=0.0) best_quantum = max((r.score for r in results if r.backend in quantum_like), default=0.0) advantage = best_quantum > best_classical + 1e-6 record = { "instance": inst.name, "n_candidates": inst.C, "n_demand_cells": inst.D, "total_risk": round(inst.total_risk(), 4), "budget_p": inst.budget_p, "reference_optimal_score": round(ref, 4), "reference_proven_optimal": bool(exact), "backends": { r.backend: { "selected": r.selected, "score": round(r.score, 4), "coverage_fraction": round(r.score / max(inst.total_risk(), 1e-9), 4), "gap_vs_reference": round(ref - r.score, 4), "latency_ms": round(r.latency_ms, 2), "proven_optimal": r.proven_optimal, "note": r.note, } for r in results }, "classical_baseline_score": round(best_classical, 4), "best_quantum_like_score": round(best_quantum, 4), "quantum_advantage_observed": bool(advantage), "fallback_available": True, "explanation": _verdict(results, ref, advantage), } return record def _verdict(results, ref, advantage) -> str: cp = next((r for r in results if r.backend == "cpsat_exact"), None) lines = [] if cp and cp.proven_optimal: lines.append(f"CP-SAT proved the optimum (risk-coverage={cp.score:.3f}) " f"in {cp.latency_ms:.1f} ms -- this is the deployable recommendation.") matched = [r.backend for r in results if r.backend != "cpsat_exact" and abs(r.score - ref) < 1e-6] if matched: lines.append("Backends matching the optimum: " + ", ".join(matched) + ".") if not advantage: lines.append("No quantum/quantum-inspired advantage observed at this scale " "(same score, higher latency than CP-SAT). Quantum stays a " "benchmark backend, not a dependency.") else: lines.append("A quantum-like backend EXCEEDED the best classical score -- " "re-verify the CP-SAT model and instance before claiming advantage.") return " ".join(lines)