""" PHOENIX instance: Sicily ground-sensor placement (Alessandria della Rocca + expansion). THIS IS A PHOENIX-ONLY INSTANCE. It supplies its own geometry, risk layer, and constraints to the neutral geo-optimize library. It shares no data with any other product instance. PROVENANCE / HONESTY BANNER --------------------------- The candidate sites, elevations, and risk weights below are SYNTHETIC but plausible — a runnable stand-in until Gaetano supplies the real artifacts: - real candidate sites (towers, ridgelines, accessible high points, existing Pi5/Hailo nodes) as GeoJSON points with elevation + mast height, - a real risk layer (EFFIS fire-danger / fuel load / ignition history / WUI asset exposure over the Sicilian AOI) on the demand grid, - real LoRa reachability + road-access flags per candidate. Swap `build_phoenix_sicily(...data_source=...)` to load those and the model is unchanged — that is the whole point of the shared formulation. Coordinates are anchored on Alessandria della Rocca, Sicily (~37.5736 N, 13.4533 E). A local equirectangular ENU frame (metres) is used for geometry; lat/lon is reconstructed for the output record. """ from __future__ import annotations import sys, os import numpy as np sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from geo_optimize import PlacementInstance ADR_LAT, ADR_LON = 37.5736, 13.4533 # Alessandria della Rocca M_PER_DEG_LAT = 111_320.0 def _enu_to_latlon(xy): lat = ADR_LAT + xy[:, 1] / M_PER_DEG_LAT lon = ADR_LON + xy[:, 0] / (M_PER_DEG_LAT * np.cos(np.radians(ADR_LAT))) return np.stack([lat, lon], axis=1) def _synthetic_terrain(gx, gy, rng): """A few ridges so line-of-sight actually matters (synthetic).""" z = np.zeros_like(gx, dtype=float) for _ in range(4): cx, cy = rng.uniform(-1, 1, 2) * 5000 ang = rng.uniform(0, np.pi) width = rng.uniform(800, 1800) d = np.abs((gx - cx) * np.cos(ang) + (gy - cy) * np.sin(ang)) z += rng.uniform(300, 700) * np.exp(-(d / width) ** 2) z += 200 + 0.01 * gy # gentle regional slope return z def build_phoenix_sicily( aoi_km: float = 12.0, demand_grid: int = 16, # 16x16 = 256 demand cells n_candidates: int = 26, budget_p: int = 10, detect_range_m: float = 5000.0, mast_m: float = 15.0, lora_range_m: float = 6500.0, qpu_scale: bool = False, # downsize so SB/QAOA/real-QPU stay small seed: int = 7, ) -> PlacementInstance: if qpu_scale: demand_grid, n_candidates, budget_p, aoi_km = 8, 12, 5, 8.0 rng = np.random.default_rng(seed) half = aoi_km * 1000 / 2 # demand grid (WUI / wildland risk cells) lin = np.linspace(-half, half, demand_grid) DX, DY = np.meshgrid(lin, lin) demand = np.stack([DX.ravel(), DY.ravel()], axis=1) z_demand = _synthetic_terrain(demand[:, 0], demand[:, 1], np.random.default_rng(seed + 1)) # synthetic risk: higher on dry south-facing mid-elevation slopes near WUI wui = np.exp(-((demand[:, 0] + 1500) ** 2 + (demand[:, 1] - 1000) ** 2) / (2 * 4000 ** 2)) fuel = np.clip((z_demand - z_demand.min()) / (np.ptp(z_demand) + 1e-9), 0, 1) ignition = rng.uniform(0.3, 1.0, demand.shape[0]) weights = (0.5 * wui + 0.3 * fuel + 0.2 * ignition) weights = weights / weights.sum() * 100.0 # normalise: total risk = 100 # candidate sites: prefer high ground (ridge tops) + the town gateway cand_pool = demand.copy() z_pool = z_demand.copy() order = np.argsort(-z_pool) # high ground first pick = order[: max(n_candidates * 3, n_candidates)] chosen = rng.choice(pick, size=n_candidates - 1, replace=False) candidates = np.vstack([cand_pool[chosen], np.array([[0.0, 0.0]])]) # last = town gateway z_cand = np.concatenate([z_pool[chosen], [z_demand[len(z_demand) // 2]]]) gateway_idx = n_candidates - 1 # coverage = within range AND coarse terrain line-of-sight C, D = n_candidates, demand.shape[0] cover = np.zeros((C, D), bool) for c in range(C): sc = candidates[c] sz = z_cand[c] + mast_m dvec = demand - sc dist = np.hypot(dvec[:, 0], dvec[:, 1]) in_range = dist <= detect_range_m for d in np.where(in_range)[0]: cover[c, d] = _los(sc, sz, demand[d], z_demand[d], demand, z_demand) # cost proxy: maintenance ~ distance from gateway (road access proxy) cost = 1.0 + np.hypot(candidates[:, 0], candidates[:, 1]) / 1000.0 # eligibility: LoRa-reachable to gateway (mesh connectivity proxy) dist_gw = np.hypot(candidates[:, 0] - candidates[gateway_idx, 0], candidates[:, 1] - candidates[gateway_idx, 1]) allowed = dist_gw <= lora_range_m allowed[gateway_idx] = True inst = PlacementInstance( candidates=candidates, demand=demand, weights=weights, cover=cover, cost=cost, budget_p=budget_p, cost_cap=None, allowed=allowed, name="phoenix_sicily_adr" + ("_qpu" if qpu_scale else ""), meta={ "product": "PHOENIX", "aoi": "Alessandria della Rocca, Sicily", "anchor_latlon": [ADR_LAT, ADR_LON], "candidates_latlon": _enu_to_latlon(candidates).tolist(), "gateway_index": gateway_idx, "detect_range_m": detect_range_m, "lora_range_m": lora_range_m, "data_provenance": "SYNTHETIC illustrative — replace with Gaetano's real " "sites + EFFIS/fuel risk layer + LoRa/road flags", "units": "local ENU metres; risk normalised to total=100", }, ) return inst def _los(src_xy, src_z, dst_xy, dst_z, grid_xy, grid_z, n_samples=6): """Coarse terrain line-of-sight: sample the straight ray and require the sight line to clear the nearest-terrain height at each sample.""" for t in np.linspace(0.15, 0.85, n_samples): p = src_xy + t * (dst_xy - src_xy) sight_z = src_z + t * (dst_z - src_z) k = np.argmin((grid_xy[:, 0] - p[0]) ** 2 + (grid_xy[:, 1] - p[1]) ** 2) if grid_z[k] > sight_z + 5.0: # terrain blocks the line return False return True if __name__ == "__main__": inst = build_phoenix_sicily() print(f"{inst.name}: C={inst.C} candidates, D={inst.D} demand cells, " f"budget={inst.budget_p}, total_risk={inst.total_risk():.1f}") print("mean cover degree per site:", inst.cover.sum(axis=1).mean().round(1)) print("demands covered by >=1 site:", (inst.cover.any(axis=0)).mean().round(3))