back to working after bend-aware straights

.gitignore

@@ -10,3 +10,5 @@ wheels/
 .venv
 
 .hypothesis
+*.png
+

README.md

@@ -67,16 +67,9 @@ if results["net1"].is_valid:
 
 ## Usage Examples
 
-Check the `examples/` directory for ready-to-run scripts demonstrating core features:
+For detailed visual demonstrations and architectural deep-dives, see the **[Examples README](examples/README.md)**.
 
-*   **`examples/01_simple_route.py`**: Basic single-net routing with visualization. Generates `01_simple_route.png`.
-*   **`examples/02_congestion_resolution.py`**: Multi-net routing resolving bottlenecks using Negotiated Congestion. Generates `02_congestion_resolution.png`.
-*   **`examples/03_locked_paths.py`**: Incremental workflow using `lock_net()` to route around previously fixed paths. Generates `03_locked_paths.png`.
-*   **`examples/04_sbends_and_radii.py`**: Complex paths using parametric S-bends and multiple bend radii. Generates `04_sbends_and_radii.png`.
-*   **`examples/05_orientation_stress.py`**: Stress test for various port orientation combinations (U-turns, opposite directions). Generates `05_orientation_stress.png`.
-*   **`examples/06_bend_collision_models.py`**: Comparison of different collision models for bends (Arc vs. BBox vs. Clipped BBox). Generates `06_bend_collision_models.png`.
-
-Run an example:
+Check the `examples/` directory for ready-to-run scripts. To run an example:
 ```bash
 python3 examples/01_simple_route.py
 ```

docs/plans/cost_and_collision_engine.md

@@ -0,0 +1,42 @@
+# Cost and Collision Engine Spec
+
+This document describes the methods for ensuring "analytic correctness" while maintaining a computationally efficient cost function.
+
+## 1. Analytic Correctness
+The router balances speed and verification using a two-tier approach.
+
+### 1.1. R-Tree Geometry Engine
+The router uses an **R-Tree of Polygons** for all geometric queries.
+*   **Move Validation:** Every "Move" proposed by the search is checked for intersection against the R-Tree.
+*   **Pre-dilation for Clearance:** All obstacles and paths use a global clearance $C$. At the start of the session, all user-provided obstacles are **pre-dilated by $(W_{max} + C)/2$** for initial broad pruning. However, individual path dilation for intersection tests uses $(W_i + C)/2$ for a specific net's width $W_i$.
+*   **Safety Zone:** To prevent immediate collision flags for ports placed on or near obstacle boundaries, the router ignores collisions within a radius of **2nm** of start and end ports.
+*   **Single Layer:** All routing and collision detection occur on a single layer.
+
+### 1.2. Cost Calculation (Soft Constraint)
+The "Danger Cost" $g_{proximity}(n)$ is a function of the distance $d$ to the nearest obstacle:
+$$g_{proximity}(n) = \begin{cases} \infty & \text{if } d < (W_i + C)/2 \\ \frac{k}{d^2} & \text{if } (W_i + C)/2 \le d < \text{Safety Threshold} \\ 0 & \text{if } d \ge \text{Safety Threshold} \end{cases}$$
+
+To optimize A* search, a **Static Danger Map (Precomputed Grid)** is used for the heuristic.
+*   **Grid Resolution:** Default **1000nm (1µm)**.
+*   **Static Nature:** The grid only accounts for fixed obstacles. It is computed once at the start of the session and is **not re-computed** during the Negotiated Congestion loop.
+*   **Efficiency:** For a 20x20mm layout, this results in a 20k x 20k matrix.
+    *   **Memory:** Using a `uint8` or `float16` representation, this consumes ~400-800MB (Default < 2GB). For extremely high resolution or larger areas, the system supports up to **20GB** allocation.
+*   **Precision:** Strict intersection checks still use the R-Tree for "analytic correctness."
+
+## 2. Collision Detection Implementation
+The system relies on `shapely` for geometry and `rtree` for spatial indexing.
+1.  **Arc Resolution:** Arcs are represented as polygons approximated by segments with a maximum deviation (sagitta).
+2.  **Intersection Test:** A "Move" is valid only if its geometry does not intersect any obstacle in the R-Tree.
+3.  **Self-Intersection:** Paths from the same net must not intersect themselves.
+4.  **No Crossings:** Strictly 2D; no crossings (vias or bridges) are supported.
+
+## 3. Negotiated Congestion (Path R-Tree)
+To handle multiple nets, the router maintains a separate R-Tree containing the dilated geometries ($C/2$) of all currently routed paths.
+*   **Congestion Cost:** $P \times (\text{Overlaps in Path R-Tree})$.
+*   **Failure Policy:** If no collision-free path is found after the max iterations, the router returns the **"least-bad" (lowest cost)** path. These paths MUST be explicitly flagged as invalid (e.g., via an `is_valid=False` attribute or a non-zero `collision_count`) so the user can identify and manually fix the failure.
+
+## 4. Handling Global Clearances
+Clearances are global. Both obstacles and paths are pre-dilated once. This ensures that any two objects maintain at least $C$ distance if their dilated versions do not intersect.
+
+## 5. Locked Paths
+The router supports **Locked Paths**—existing geometries inserted into the static Obstacle R-Tree, ensuring they are never modified or rerouted.

docs/plans/geometric_representation.md

@@ -0,0 +1,37 @@
+# Geometric Representation and Data Structures
+
+This document defines the core data structures for representing the routing problem.
+
+## 1. Port Definitions (Connectivity)
+Routing requests are defined as a mapping of source ports to destination ports.
+
+*   **Port Structure:** `(x, y, orientation)`
+    *   `x, y`: Coordinates snapped to **1nm** grid.
+    *   `orientation`: Strictly $\{0, 90, 180, 270\}$ degrees.
+*   **Netlist:** A dictionary or list of tuples: `{(start_port): (end_port)}`. Each entry represents a single-layer path that must be routed.
+*   **Heterogeneous Widths:** The router supports different widths $W_i$ per net. Dilation for a specific net $i$ is calculated as $(W_i + C)/2$ to maintain the global clearance $C$.
+
+## 2. Component Library (Move Generator)
+The router uses a discrete set of components to expand states in the A* search.
+
+### 2.1. Straight Waveguides
+*   **A* Expansion:** Generates "Straight" moves of varying lengths (e.g., $1\mu m, 5\mu m, 25\mu m$).
+*   **Snap-to-Target:** If the destination port $T$ is directly ahead and the current state's orientation matches $T$'s entry orientation, a special move is generated to close the gap exactly.
+
+### 2.2. Fixed 90° Bends (PDK Cells)
+*   **Parameters:** `radius`, `width`.
+*   **A* Expansion:** A discrete move that changes orientation by $\pm 90^\circ$ and shifts the coordinate by the radius.
+*   **Grid Alignment:** If the bend radius $R$ is not a multiple of the search grid (default $1\mu m$), the resulting state is **snapped to the nearest grid point**, and a warning is issued to the user.
+
+### 2.3. Parametric S-Bends (Compact)
+*   **Parameters:** `radius`, `width`.
+*   **A* Expansion:** Used ONLY for lateral offsets $O < 2R$.
+*   **Large Offsets ($O \ge 2R$):** The router does not use a single S-bend move for large offsets. Instead, the A* search naturally finds the optimal path by combining two 90° bends and a straight segment. This ensures maximum flexibility in obstacle avoidance for large shifts.
+
+## 3. Obstacle Representation
+Obstacles are provided as raw polygons on a single layer.
+
+*   **Pre-processing:** All input polygons are inserted into an **R-Tree**.
+*   **Buffer/Dilation:** Obstacles are pre-dilated by $(W_{max} + Clearance)/2$ for initial pruning, but final collision tests use the net-specific width $W_i$.
+*   **No Multi-layer:** The router assumes all obstacles and paths share the same plane.
+*   **Safety Zone:** Ignore collisions within **2nm** of start and end ports for robustness.

docs/plans/high_level_notes.md

@@ -0,0 +1,13 @@
+# High-Level Notes: Auto-Routing for Integrated Circuits
+
+We're implementing auto-routing for photonic and RF integrated circuits. The problem space has the following features:
+
+    - **Single Layer:** All paths are on a single layer; multi-layer routing is not supported.
+    - **No Crossings:** Crossings are to be avoided and are not supported by the router's automatic placement. The user must manually handle any required crossings.
+    - **Large Bend Radii:** Bends use large radii ($R$), and are usually pre-generated (e.g. 90-degree cells).
+    - **Proximity Sensitivity:** Paths are sensitive to proximity to other paths and obstacles (crosstalk/coupling).
+    - **Manhattan Preference:** Manhattan pathing is sufficient for most cases; any-angle is rare.
+    - **S-Bends:** S-bends are necessary for small lateral offsets ($O < 2R$).
+    - **Hybrid Search:** A* state-lattice search is used for discrete component placement.
+    - **Dilation:** A $Clearance/2$ dilation is applied to all obstacles and paths for efficient collision avoidance.
+    - **Negotiated Congestion:** A multi-net "PathFinder" loop iteratively reroutes nets to resolve congestion.

docs/plans/implementation_plan.md

@@ -0,0 +1,64 @@
+# Implementation Plan
+
+This plan outlines the step-by-step implementation of the `inire` auto-router. For detailed test cases, refer to [Testing Plan](./testing_plan.md).
+
+## Phase 1: Core Geometry & Move Generation
+**Goal:** Implement Ports, Polygons, and Component Library with high geometric fidelity.
+1.  **Project Setup:** Initialize `inire/` structure and `pytest` configuration. Include `hypothesis` for property-based testing.
+2.  **`geometry.primitives`:**
+    *   `Port` with **1nm** snapping.
+    *   Basic 2D transformations (rotate, translate).
+    *   **Property-Based Tests:** Verify transform invariants (e.g., $90^\circ$ rotation cycles).
+3.  **`geometry.components`:**
+    *   `Straight`, `Bend90`, `SBend`.
+    *   **Search Grid Snapping:** Implement 1µm snapping for expanded ports.
+    *   **Small S-Bends ($O < 2R$):** Logic for parametric generation.
+    *   **Edge Cases:** Handle $O=2R$ and $L < 1\mu m$.
+4.  **Tests:** 
+    *   Verify geometric correctness (refer to Testing Plan Section 1).
+    *   Unit tests for `Port` snapping and component transformations.
+
+## Phase 2: Collision Engine & Cost
+**Goal:** Build the R-Tree wrapper and the analytic cost function.
+1.  **`geometry.collision`:** Implement `CollisionEngine`.
+    *   **Pre-dilation:** Obstacles/Paths dilated by $Clearance/2$.
+    *   **Safety Zone:** Ignore collisions within **2nm** of start/end ports.
+2.  **`router.danger_map`:**
+    *   Implement **1µm** pre-computed proximity grid.
+    *   Optimize for design sizes up to **20x20mm** (< 2GB memory).
+3.  **`router.cost`:** Implement `CostEvaluator`.
+    *   Bend cost: $10 \times (\text{Manhattan distance between ports})$.
+    *   Integrate R-Tree for strict checks and Danger Map for heuristic.
+4.  **Tests:** 
+    *   Verify collision detection with simple overlapping shapes (Testing Plan Section 2.1).
+    *   Verify Danger Map accuracy and memory footprint (Testing Plan Section 2.2).
+    *   **Post-Route Validator:** Implement the independent `validate_path` utility.
+
+## Phase 3: Single-Net A* Search
+**Goal:** Route a single net from A to B with 1nm precision.
+1.  **`router.astar`:** Implement the priority queue loop.
+    *   State representation: `(x_µm, y_µm, theta)`.
+    *   Move expansion loop with 1µm grid.
+    *   **Natural S-Bends:** Ensure search can find $O \ge 2R$ shifts by combining moves.
+    *   **Look-ahead Snapping:** Actively bridge to the 1nm target when in the capture radius (10µm).
+2.  **Heuristic:** Manhattan distance $h(n)$ + orientation penalty + Danger Map lookup.
+3.  **Tests:** 
+    *   Solve simple maze problems and verify path optimality (Testing Plan Section 3).
+    *   Verify snap-to-target precision at 1nm resolution.
+    *   **Determinism:** Verify same seed = same path.
+
+## Phase 4: Multi-Net PathFinder
+**Goal:** Implement the "Negotiated Congestion" loop for multiple nets.
+1.  **`router.pathfinder`:**
+    *   Sequential routing -> Identify congestion -> Inflate cost -> Reroute.
+    *   **R-Tree Congestion:** Store dilated path geometries.
+2.  **Explicit Results:** Return `RoutingResult` objects with `is_valid` and `collisions` metadata.
+3.  **Tests:** 
+    *   Full multi-net benchmarks (Testing Plan Section 4).
+    *   Verify rerouting behavior in crowded environments.
+
+## Phase 5: Visualization, Benchmarking & Fuzzing
+1.  **`utils.visualization`:** Plot paths using `matplotlib`. Highlight collisions in red.
+2.  **Benchmarks:** Stress test with 50+ nets. Verify performance and node limits (Testing Plan Section 5).
+3.  **Fuzzing:** Run A* on randomized layouts to ensure stability.
+4.  **Final Validation:** Ensure all `is_valid=True` results pass the independent `validate_path` check.

docs/plans/package_structure.md

@@ -0,0 +1,57 @@
+# Python Package Structure
+
+This document outlines the directory structure and module organization for the `inire` auto-router package.
+
+## 1. Directory Layout
+
+```
+inire/
+├── __init__.py             # Exposes the main `Router` class and key types
+├── geometry/               # Core geometric primitives and operations
+│   ├── __init__.py
+│   ├── primitives.py       # Point, Port, Polygon, Arc classes
+│   ├── collision.py        # R-Tree wrapper and intersection logic
+│   └── components.py       # Move generators (Straight, Bend90, SBend)
+├── router/                 # Search algorithms and pathfinding
+│   ├── __init__.py
+│   ├── astar.py            # Hybrid State-Lattice A* implementation
+│   ├── graph.py            # Node, Edge, and Graph data structures
+│   ├── cost.py             # Cost functions (length, bend, proximity)
+│   ├── danger_map.py       # Pre-computed grid for heuristic proximity costs
+│   └── pathfinder.py       # Multi-net "Negotiated Congestion" manager
+├── utils/                  # Utility functions
+│   ├── __init__.py
+│   └── visualization.py    # Plotting tools for debug/heatmaps (matplotlib/klayout)
+└── tests/                  # Unit and integration tests
+    ├── __init__.py
+    ├── conftest.py         # Pytest fixtures (common shapes, PDK cells)
+    ├── test_primitives.py  # Tests for Port and coordinate transforms
+    ├── test_components.py  # Tests for Straight, Bend90, SBend generation
+    ├── test_collision.py   # Tests for R-Tree and dilation logic
+    ├── test_cost.py        # Tests for Danger Map and cost evaluation
+    ├── test_astar.py       # Tests for single-net routing (mazes, snapping)
+    └── test_pathfinder.py  # Multi-net "Negotiated Congestion" benchmarks
+```
+
+## 2. Module Responsibilities
+
+### `inire.geometry`
+*   **`primitives.py`**: Defines the `Port` named tuple `(x, y, theta)` and helper functions for coordinate transforms.
+*   **`collision.py`**: Wraps the `rtree` or `shapely` library. Handles the "Analytic Correctness" checks (exact polygon distance).
+*   **`components.py`**: Logic to generate "Moves" from a start port. E.g., `SBend.generate(start_port, offset, radius)` returns a list of polygons and the end port. Handles $O > 2R$ logic.
+
+### `inire.router`
+*   **`astar.py`**: The heavy lifter. Maintains the `OpenSet` (priority queue) and `ClosedSet`. Implements the "Snap-to-Target" logic.
+*   **`cost.py`**: compute $f(n) = g(n) + h(n)$. encapsulates the "Danger Map" and Path R-Tree lookups.
+*   **`danger_map.py`**: Manages the pre-computed proximity grid used for $O(1)$ heuristic calculations.
+*   **`pathfinder.py`**: Orchestrates the multi-net loop. Tracks the Path R-Tree for negotiated congestion and triggers reroutes.
+
+### `inire.tests`
+*   **Structure:** Tests are co-located within the package for ease of access.
+*   **Fixtures:** `conftest.py` will provide standard PDK cells (e.g., a $10\mu m$ radius bend) to avoid repetition in test cases.
+
+## 3. Dependencies
+*   `numpy`: Vector math.
+*   `shapely`: Polygon geometry and intersection.
+*   `rtree`: Spatial indexing.
+*   `networkx` (Optional): Not used for core search to ensure performance.

docs/plans/routing_architecture_decision.md

@@ -0,0 +1,67 @@
+# Architecture Decision Record: Auto-Routing for Photonic & RF ICs
+
+## 1. Problem Context
+Photonic and RF routing differ significantly from digital VLSI due to physical constraints:
+* **Geometric Rigidity:** 90° bends are pre-rendered PDK cells with fixed bounding boxes.
+* **Parametric Flexibility:** S-bends must be generated on-the-fly to handle arbitrary offsets, provided they maintain a constant radius $R$.
+* **Signal Integrity:** High sensitivity to proximity (coupling/crosstalk) and a strong preference for single-layer, non-crossing paths.
+* **Manual Intervention:** The router is strictly 2D and avoids all other geometries on the same layer. No crossings (e.g. vias or bridges) are supported by the automatic routing engine. The user must manually handle any required crossings by placing components (e.g. crossing cells) and splitting the net list accordingly. This simplifies the router's task to 2D obstacle avoidance and spacing optimization.
+
+---
+
+## 2. Candidate Algorithms & Tradeoffs
+
+### 2.1. Rubber-Band (Topological) Routing
+This approach treats paths as elastic bands stretched around obstacles, later "inflating" them to have width and curvature.
+
+| Feature | Analysis |
+| :--- | :--- |
+| **Strengths** | Excellent at "River Routing" and maintaining minimum clearances. Inherently avoids crossings. |
+| **Downsides** | **The Inflation Gap:** A valid thin-line topology may be physically un-routable once the large radius $R$ is applied. It struggles to integrate rigid, pre-rendered 90° blocks into a continuous elastic model. |
+| **Future Potential** | High, if a "Post-Processing" engine can reliably snap elastic curves to discrete PDK cells without breaking connectivity. |
+
+### 2.2. Voronoi-Based (Medial Axis) Routing
+Uses a Voronoi diagram to find paths that are maximally distant from all obstacles.
+
+| Feature | Analysis |
+| :--- | :--- |
+| **Strengths** | Theoretically optimal for "Safety" and crosstalk reduction. Guaranteed maximum clearance. |
+| **Downsides** | **Manhattan Incompatibility:** Voronoi edges are any-angle and often jagged. Mapping these to a Manhattan-heavy PDK (90° bends) requires a lossy "snapping" phase that often violates the very safety the algorithm intended to provide. |
+| **Future Potential** | Useful as a "Channel Finder" to guide a more rigid router, but unsuitable as a standalone geometric solver. |
+
+### 2.3. Integer Linear Programming (ILP)
+Formulates routing as a massive optimization problem where a solver picks the best path from a pool of candidates.
+
+| Feature | Analysis |
+| :--- | :--- |
+| **Strengths** | Can find the mathematically "best" layout (e.g., minimum total length or total bends) across all nets simultaneously. |
+| **Downsides** | **Candidate Explosion:** Because S-bends are generated on-the-fly, the number of possible candidate shapes is infinite. To use ILP, one must "discretize" the search space, which may miss the one specific geometry needed for a tight fit. |
+| **Future Potential** | Effective for small, high-congestion "Switchbox" areas where all possible geometries can be pre-tabulated. |
+
+---
+
+## 3. Selected Approach: Hybrid State-Lattice A*
+
+### 3.1. Rationale
+The **State-Lattice** variant of the A* algorithm is selected as the primary routing engine. Unlike standard A* which moves between grid cells, State-Lattice moves between **states** defined as $(x, y, \theta)$.
+
+1.  **Native PDK Integration:** The router treats the pre-rendered 90° bend cell as a discrete "Move" in the search tree. The algorithm only considers placing a bend if the cell's bounding box is clear of obstacles.
+2.  **Parametric S-Bends:** When the search needs to bridge a lateral gap, it triggers a "Procedural Move." It calculates a fixed-radius S-bend on-the-fly. If the resulting arc is valid and collision-free, it is added as an edge in the search graph.
+3.  **Predictable Costing:** It allows for a sophisticated cost function $f(n) = g(n) + h(n)$ where:
+    * $g(n)$ penalizes path length and proximity to obstacles (using a distance-transform field).
+    * $h(n)$ (the heuristic) guides the search toward the destination while favoring Manhattan alignments.
+
+### 3.2. Implementation Strategy
+* **Step 1: Distance Transform.** Pre-calculate a "Danger Map" of the layout. Cells close to obstacles have a high cost; cells far away have low cost. This handles the **Proximity Sensitivity** constraint.
+* **Step 2: State Expansion.** From the current point, explore:
+    * `Straight(length)`
+    * `PDK_Bend_90(direction)`
+    * `S_Bend(target_offset, R)`
+* **Step 3: Rip-up and Reroute.** To handle the sequential nature of A*, implement a "Negotiated Congestion" scheme (PathFinder algorithm) where nets "pay" more to occupy areas that other nets also want.
+
+---
+
+## 4. Summary of Tradeoffs for Future Review
+* **Why not pure Manhattan?** Photonic/RF requirements for large $R$ and S-bends make standard grid-based maze routers obsolete.
+* **Why not any-angle?** While any-angle is possible, the PDK's reliance on pre-rendered 90° blocks makes a lattice-based approach more manufacturing-stable.
+* **Risk:** The primary risk is **Search Time**. As the library of moves grows (more S-bend options), the branching factor increases. This must be managed with aggressive pruning and spatial indexing (e.g., R-trees).

docs/plans/routing_search_spec.md

@@ -0,0 +1,66 @@
+# Routing Search Specification
+
+This document details the Hybrid State-Lattice A* implementation and the "PathFinder" (Negotiated Congestion) algorithm for multi-net routing.
+
+## 1. Hybrid State-Lattice A* Search
+The core router operates on states $S = (x, y, \theta)$, where $\theta \in \{0, 90, 180, 270\}$.
+
+### 1.1. State Representation & Snapping
+To ensure search stability and hash-map performance:
+*   **Intermediate Ports:** Every state expanded by A* is snapped to a search grid.
+*   **Search Grid:** Default snap is **1000nm (1µm)**.
+*   **Final Alignment:** The "Snap-to-Target" logic bridges the gap from the coarse search grid to the final **1nm** port coordinates.
+*   **Max Design Size:** Guidelines for memory/performance assume up to **20mm x 20mm** routing area.
+
+### 1.2. Move Expansion Logic
+From state $S_n$, the following moves are evaluated:
+1.  **Tiered Straight Steps:** Expand by a set of distances $L \in \{1\mu m, 5\mu m, 25\mu m\}$.
+2.  **Snap-to-Target:** A "last-inch" look-ahead move. If the target $T$ is within a **Capture Radius (Default: 10µm)** and a straight segment or single bend can reach it, a special move is generated to close the gap exactly at 1nm resolution.
+3.  **90° Bend:** Try clockwise/counter-clockwise turns using fixed PDK cells.
+4.  **Small-Offset S-Bend:**
+    *   **Only for $O < 2R$:** Parametric S-bend (two tangent arcs).
+    *   **$O \ge 2R$:** Search naturally finds these by combining 90° bends and straight segments.
+
+### 1.3. Cost Function $f(n) = g(n) + h(n)$
+The search uses a flexible, component-based cost model.
+
+*   **$g(n)$ (Actual Cost):** $\sum \text{ComponentCost}_i + \text{ProximityCost} + \text{CongestionCost}$
+    *   **Straight Segment:** $L \times C_{unit\_length}$.
+    *   **90° Bend:** $10 \times (\text{Manhattan distance between ports})$.
+    *   **S-Bend:** $f(O, R)$.
+    *   **Proximity Cost:** $k/d^2$ penalty (strict checks use R-Tree).
+    *   **Congestion Cost:** $P \times (\text{Overlaps in Path R-Tree})$.
+*   **$h(n)$ (Heuristic):**
+    *   Manhattan distance $L_1$ to the target.
+    *   Orientation Penalty: High cost if the state's orientation doesn't match the target's entry orientation.
+    *   **Greedy Weighting:** The A* search uses a weighted heuristic (e.g., $1.1 \times h(n)$) to prioritize search speed over strict path optimality.
+    *   **Danger Map Heuristic:** Fast lookups from the **1µm** pre-computed proximity grid.
+
+## 2. Multi-Net "PathFinder" Strategy (Negotiated Congestion)
+1.  **Iteration:** Identify "Congestion Areas" using Path R-Tree intersections.
+2.  **Inflation:** Increase penalty multiplier $P$ for congested areas.
+3.  **Repeat:** Continue until no overlaps exist or the max iteration count is reached (Default: **20 iterations**).
+
+### 2.1. Convergence & Result Policy
+*   **Least Bad Attempt:** If no 100% collision-free solution is found, return the result with the lowest total cost (including overlaps).
+*   **Explicit Reporting:** Results MUST include a `RoutingResult` object containing:
+    *   `path_geometry`: The actual polygon sequence.
+    *   `is_valid`: Boolean (True only if no collisions).
+    *   `collisions`: A count or list of detected overlap points/polygons.
+*   **Visualization:** Overlapping regions are highlighted (e.g., in red) in the heatmaps.
+
+## 3. Search Limits & Scaling
+*   **Node Limit:** A* search is capped at **50,000 nodes** per net per iteration.
+*   **Dynamic Timeout:** Session-level timeout based on problem size:
+    *   `Timeout = max(2s, 0.05s * num_nets * num_iterations)`.
+    *   *Example:* A 100-net problem over 20 iterations times out at **100 seconds**.
+
+## 4. Determinism
+All search and rip-up operations are strictly deterministic.
+*   **Seed:** A user-provided `seed` (int) MUST be used to initialize any random number generators (e.g., if used for tie-breaking or initial net ordering).
+*   **Tie-Breaking:** If two nodes have the same $f(n)$, a consistent tie-breaking rule (e.g., based on node insertion order or state hash) must be used.
+
+## 5. Optimizations
+*   **A* Pruning:** Head toward the target and prune suboptimal orientations.
+*   **Safety Zones:** Ignore collisions within **2nm** of start/end ports.
+*   **Spatial Indexing:** R-Tree queries are limited to the bounding box of the current move.

docs/plans/testing_plan.md

@@ -0,0 +1,138 @@
+# Testing Plan (Refined)
+
+This document defines the comprehensive testing strategy for the `inire` auto-router, ensuring analytic correctness and performance.
+
+## 1. Unit Tests: Geometry & Components (`inire/geometry`)
+
+### 1.1. Primitives (`primitives.py`)
+*   **Port Snapping:** Verify that `Port(x, y, orientation)` snaps `x` and `y` to the nearest **1nm** grid point.
+*   **Coordinate Transforms:**
+    *   Translate a port by `(dx, dy)`.
+    *   Rotate a port by `90`, `180`, `270` degrees around an origin.
+    *   Verify orientation wrapping (e.g., $270^\circ + 90^\circ \to 0^\circ$).
+*   **Polygon Creation:**
+    *   Generate a polygon from a list of points.
+    *   Verify bounding box calculation.
+*   **Property-Based Testing (Hypothesis):**
+    *   Verify that any `Port` after transformation remains snapped to the **1nm** grid.
+    *   Verify that $Rotate(Rotate(Port, 90), -90)$ returns the original `Port` (up to snapping).
+
+### 1.2. Component Generation (`components.py`)
+*   **Straight Moves:**
+    *   Generate a `Straight(length=10.0, width=2.0)` segment.
+    *   Verify the end port is exactly `length` away in the correct orientation.
+    *   Verify the resulting polygon's dimensions.
+    *   **Edge Case:** $L < 1\mu m$ (below search grid). Verify it still generates a valid 1nm segment.
+*   **Bend90:**
+    *   Generate a `Bend90(radius=10.0, width=2.0, direction='CW')`.
+    *   Verify the end port's orientation is changed by $-90^\circ$.
+    *   Verify the end port's position is exactly `(R, -R)` relative to the start (for $0^\circ \to 270^\circ$).
+    *   **Grid Snapping:** Verify that if the bend radius results in a non-1µm aligned port, it is snapped to the nearest **1µm** search grid point (with a warning).
+*   **SBend (Small Offset $O < 2R$):**
+    *   Generate an `SBend(offset=5.0, radius=10.0, width=2.0)`.
+    *   Verify the total length matches the analytical $L = 2\sqrt{O(2R - O/4)}$ (or equivalent arc-based formula).
+    *   Verify the tangent continuity at the junction of the two arcs.
+    *   **Edge Case:** $O = 2R$. Verify it either generates two 90-degree bends or fails gracefully with a clear error.
+    *   Verify it fails/warns if $O > 2R$.
+
+### 1.3. Geometric Fidelity
+*   **Arc Resolution (Sagitta):**
+    *   Verify that `Bend90` and `SBend` polygons are approximated by segments such that the maximum deviation (sagitta) is within a user-defined tolerance (e.g., 10nm).
+    *   Test with varying radii to ensure segment count scales appropriately.
+
+## 2. Unit Tests: Collision & Cost (`inire/geometry/collision` & `router/cost`)
+
+### 2.1. Collision Engine
+*   **Pre-dilation Logic:**
+    *   Verify that an obstacle (polygon) is correctly dilated by $(W_{max} + C)/2$.
+    *   **Heterogeneous Widths:** Verify that a path for Net A (width $W_1$) is dilated by $(W_1 + C)/2$, while Net B (width $W_2$) uses $(W_2 + C)/2$.
+*   **Locked Paths:**
+    *   Insert an existing path geometry into the "Static Obstacle" R-Tree.
+    *   Verify that the router treats it as an unmovable obstacle and avoids it.
+*   **R-Tree Queries:**
+    *   Test intersection detection between two overlapping polygons.
+    *   Test non-intersection between adjacent but non-overlapping polygons (exactly $C$ distance apart).
+*   **Safety Zone (2nm):**
+    *   Create a port exactly on the edge of an obstacle.
+    *   Verify that a "Move" starting from this port is NOT flagged for collision if the intersection occurs within **2nm** of the port.
+*   **Self-Intersection:**
+    *   Verify that a path consisting of multiple segments is flagged if it loops back on itself.
+
+### 2.2. Danger Map & Cost Evaluator
+*   **Danger Map Generation:**
+    *   Initialize a map for a $100\mu m \times 100\mu m$ area with a single obstacle.
+    *   Verify the cost $g_{proximity}$ matches $k/d^2$ for cells near the obstacle.
+    *   Verify cost is $0$ for cells beyond the **Safety Threshold**.
+*   **Memory Check:**
+    *   Mock a $20mm \times 20mm$ grid and verify memory allocation stays within limits (e.g., `< 2GB` for standard `uint8` resolution).
+*   **Cost Calculation:**
+    *   Verify total cost $f(n)$ correctly sums length, bend penalties ($10 \times$ Manhattan), and proximity costs.
+
+### 2.3. Robustness & Limits
+*   **Design Bounds:**
+    *   Test routing at the extreme edges of the $20mm \times 20mm$ coordinate space.
+    *   Verify that moves extending outside the design bounds are correctly pruned or flagged.
+*   **Empty/Invalid Inputs:**
+    *   Test with an empty netlist.
+    *   Test with start and end ports at the exact same location.
+
+## 3. Integration Tests: Single-Net A* Search (`inire/router/astar`)
+
+### 3.1. Open Space Scenarios
+*   **Straight Line:** Route from `(0,0,0)` to `(100,0,0)`. Verify it uses only `Straight` moves.
+*   **Simple Turn:** Route from `(0,0,0)` to `(20,20,90)`. Verify it uses a `Bend90` and `Straight` segments.
+*   **Small S-Bend:** Route with an offset of $5\mu m$ and radius $10\mu m$. Verify it uses the `SBend` component.
+*   **Large S-Bend ($O \ge 2R$):** Route with an offset of $50\mu m$ and radius $10\mu m$. Verify it combines two `Bend90`s and a `Straight` segment.
+
+### 3.2. Obstacle Avoidance (The "Maze" Tests)
+*   **L-Obstacle:** Place an obstacle blocking the direct path. Verify the router goes around it.
+*   **Narrow Channel:** Create two obstacles with a gap slightly wider than $W_i + C$. Verify the router passes through.
+*   **Dead End:** Create a U-shaped obstacle. Verify the search explores alternatives and fails gracefully if no path exists.
+
+### 3.3. Snapping & Precision
+*   **Snap-to-Target Lookahead:**
+    *   Route to a target at `(100.005, 0, 0)` (not on 1µm grid).
+    *   Verify the search reaches the vicinity via the 1µm grid and the final segment bridges the **5nm** gap exactly.
+*   **Grid Alignment:**
+    *   Start from a port at `(0.5, 0.5, 0)`. Verify it snaps to the 1µm search grid correctly for the first move expansion.
+
+### 3.4. Failure Modes
+*   **Unreachable Target:** Create a target completely enclosed by obstacles. Verify the search terminates after exploring all options (or hitting the 50,000 node limit) and returns an invalid result.
+*   **Start/End Collision:** Place a port deep inside an obstacle (beyond the 2nm safety zone). Verify the router identifies the immediate collision and fails gracefully.
+
+## 4. Integration Tests: Multi-Net PathFinder (`inire/router/pathfinder`)
+
+### 4.1. Congestion Scenarios
+*   **Parallel Paths:** Route two nets that can both take straight paths. Verify no reroutes occur.
+*   **The "Cross" Test:** Two nets must cross paths in 2D.
+    *   Since crossings are illegal, verify the second net finds a detour.
+    *   Verify the `Negotiated Congestion` loop increases the cost of the shared region.
+*   **Bottleneck:** Force 3 nets through a channel that only fits 2.
+    *   Verify the router returns 2 valid paths and 1 "least bad" path (with collisions flagged).
+    *   Verify the `is_valid=False` attribute is set for the failing net.
+
+### 4.2. Determinism & Performance
+*   **Seed Consistency:** Run the same multi-net problem twice with the same seed; verify identical results (pixel-perfect).
+*   **Node Limit Enforcement:** Trigger a complex search that exceeds **50,000 nodes**. Verify it terminates and returns the best-so-far or failure.
+*   **Timeout:** Verify the session-level timeout stops the process for extremely large problems.
+
+## 5. Benchmarking & Regression
+*   **Standard Benchmark Suite:** A set of 5-10 layouts with varying net counts (1 to 50).
+*   **Metrics to Track:**
+    *   Total wire length.
+    *   Total number of bends.
+    *   Execution time per net.
+    *   Success rate (percentage of nets routed without collisions).
+    *   **Node Expansion Rate:** Nodes per second.
+    *   **Memory Usage:** Peak RSS during 20x20mm routing.
+*   **Fuzz Testing:**
+    *   Generate random obstacles and ports within a 1mm x 1mm area.
+    *   Verify that the router never crashes.
+    *   Verify that every result marked `is_valid=True` is confirmed collision-free by a high-precision (slow) check.
+
+## 6. Analytic Correctness Guarantees
+*   **Post-Route Validation:**
+    *   Implement an independent `validate_path(path, obstacles, clearance)` function using `shapely`'s most precise intersection tests.
+    *   Run this on every test result to ensure the `CollisionEngine` (which uses R-Tree for speed) hasn't missed any edge cases.
+*   **Orientation Check:**
+    *   Verify that the final port of every path matches the target orientation exactly $\{0, 90, 180, 270\}$.

examples/01_simple_route.py

@@ -26,8 +26,8 @@ def main() -> None:
     # Precompute the danger map (distance field) for heuristics
     danger_map.precompute([obstacle])
 
-    evaluator = CostEvaluator(engine, danger_map)
-    router = AStarRouter(evaluator)
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    router = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     pf = PathFinder(router, evaluator)
 
     # 2. Define Netlist

examples/02_congestion_resolution.py

@@ -16,8 +16,8 @@ def main() -> None:
     danger_map = DangerMap(bounds=bounds)
     danger_map.precompute([])
 
-    evaluator = CostEvaluator(engine, danger_map)
-    router = AStarRouter(evaluator)
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    router = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     pf = PathFinder(router, evaluator)
 
     # 2. Define Netlist

examples/03_locked_paths.py

@@ -1,3 +1,5 @@
+from shapely.geometry import Polygon
+
 from inire.geometry.collision import CollisionEngine
 from inire.geometry.primitives import Port
 from inire.router.astar import AStarRouter
@@ -8,63 +10,51 @@ from inire.utils.visualization import plot_routing_results
 
 
 def main() -> None:
-    print("Running Example 03: Locked Paths (Incremental Routing - Bus Scenario)...")
+    print("Running Example 03: Locked Paths...")
 
     # 1. Setup Environment
-    bounds = (0, 0, 120, 120)
+    bounds = (0, 0, 100, 100)
     engine = CollisionEngine(clearance=2.0)
     danger_map = DangerMap(bounds=bounds)
-    danger_map.precompute([]) # Start with empty space
+    danger_map.precompute([])
 
-    evaluator = CostEvaluator(engine, danger_map, greedy_h_weight=1.2)
-    router = AStarRouter(evaluator, node_limit=200000)
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    router = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     pf = PathFinder(router, evaluator)
 
-    # 2. Phase 1: Route a "Bus" of 3 parallel nets
-    # We give them a small jog to make the locked geometry more interesting
-    netlist_p1 = {
-        "bus_0": (Port(10, 40, 0), Port(110, 45, 0)),
-        "bus_1": (Port(10, 50, 0), Port(110, 55, 0)),
-        "bus_2": (Port(10, 60, 0), Port(110, 65, 0)),
+    # 2. Add a 'Pre-routed' net and lock it
+    # Net 'fixed' goes right through the middle
+    fixed_start = Port(10, 50, 0)
+    fixed_target = Port(90, 50, 0)
+    
+    print("Routing initial net...")
+    res_fixed = router.route(fixed_start, fixed_target, net_width=2.0)
+    
+    if res_fixed:
+        # 3. Lock this net! It now behaves like a static obstacle
+        geoms = [comp.geometry[0] for comp in res_fixed]
+        engine.add_path("locked_net", geoms)
+        engine.lock_net("locked_net")
+        print("Initial net locked as static obstacle.")
+    
+    # Update danger map to reflect the new static obstacle
+    danger_map.precompute(list(engine.static_geometries.values()))
+
+    # 4. Route a new net that must detour around the locked one
+    netlist = {
+        "detour_net": (Port(50, 10, 90), Port(50, 90, 90)),
     }
-    print("Phase 1: Routing bus (3 nets)...")
-    results_p1 = pf.route_all(netlist_p1, dict.fromkeys(netlist_p1, 2.0))
+    net_widths = {"detour_net": 2.0}
 
-    # Lock all Phase 1 nets
-    path_polys = []
-    for nid, res in results_p1.items():
-        if res.is_valid:
-            print(f"  Locking {nid}...")
-            engine.lock_net(nid)
-            path_polys.extend([p for comp in res.path for p in comp.geometry])
-        else:
-            print(f"  Warning: {nid} failed to route correctly.")
-
-    # Update danger map with the newly locked geometry
-    print("Updating DangerMap with locked paths...")
-    danger_map.precompute(path_polys)
-
-    # 3. Phase 2: Route secondary nets that must navigate around the locked bus
-    # These nets cross the bus vertically.
-    netlist_p2 = {
-        "cross_left": (Port(30, 10, 90), Port(30, 110, 90)),
-        "cross_right": (Port(80, 110, 270), Port(80, 10, 270)), # Top to bottom
-    }
-
-    print("Phase 2: Routing crossing nets around locked bus...")
-    # We use a slightly different width for variety
-    results_p2 = pf.route_all(netlist_p2, dict.fromkeys(netlist_p2, 1.5))
-
-    # 4. Check Results
-    for nid, res in results_p2.items():
-        status = "Success" if res.is_valid else "Failed"
-        print(f"  {nid:12}: {status}, collisions={res.collisions}")
+    print("Routing detour net around locked path...")
+    results = pf.route_all(netlist, net_widths)
 
     # 5. Visualize
-    all_results = {**results_p1, **results_p2}
-    all_netlists = {**netlist_p1, **netlist_p2}
-
-    fig, ax = plot_routing_results(all_results, [], bounds, netlist=all_netlists)
+    # Add the locked net back to results for display
+    from inire.router.pathfinder import RoutingResult
+    display_results = {**results, "locked_net": RoutingResult("locked_net", res_fixed or [], True, 0)}
+    
+    fig, ax = plot_routing_results(display_results, list(engine.static_geometries.values()), bounds, netlist=netlist)
     fig.savefig("examples/03_locked_paths.png")
     print("Saved plot to examples/03_locked_paths.png")
 

examples/04_sbends_and_radii.py

@@ -1,4 +1,3 @@
-
 from inire.geometry.collision import CollisionEngine
 from inire.geometry.primitives import Port
 from inire.router.astar import AStarRouter
@@ -23,15 +22,14 @@ def main() -> None:
         danger_map,
         unit_length_cost=1.0,
         greedy_h_weight=1.5,
+        bend_penalty=10.0,
+        sbend_penalty=20.0,
     )
 
-    # We want a 45 degree switchover for S-bend.
-    # Offset O = 2 * R * (1 - cos(theta))
-    # If R = 10, O = 5.86
-
     router = AStarRouter(
         evaluator,
         node_limit=50000,
+        snap_size=1.0,
         bend_radii=[10.0, 30.0],
         sbend_offsets=[5.0], # Use a simpler offset
         sbend_radii=[10.0],

examples/05_orientation_stress.py

@@ -11,41 +11,31 @@ def main() -> None:
     print("Running Example 05: Orientation Stress Test...")
 
     # 1. Setup Environment
-    # Give some breathing room (-20 to 120) for U-turns and flips (R=10)
-    bounds = (-20, -20, 120, 120)
+    bounds = (0, 0, 200, 200)
     engine = CollisionEngine(clearance=2.0)
     danger_map = DangerMap(bounds=bounds)
     danger_map.precompute([])
 
-    evaluator = CostEvaluator(engine, danger_map, greedy_h_weight=1.1)
-#    router = AStarRouter(evaluator, node_limit=100000)
-    router = AStarRouter(evaluator, node_limit=100000, bend_collision_type="clipped_bbox", bend_clip_margin=1.0)
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    router = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     pf = PathFinder(router, evaluator)
 
-    # 2. Define Netlist with various orientation challenges
+    # 2. Define Netlist: Complex orientation challenges
     netlist = {
-        # Opposite directions: requires two 90-degree bends to flip orientation
-        "opposite": (Port(10, 80, 0), Port(90, 80, 180)),
-
-        # 90-degree turn: standard L-shape
-        "turn_90": (Port(10, 60, 0), Port(40, 90, 90)),
-
-        # Output behind input: requires a full U-turn
-        "behind": (Port(80, 40, 0), Port(20, 40, 0)),
-
-        # Sharp return: output is behind and oriented towards the input
-        "return_loop": (Port(80, 20, 0), Port(40, 10, 180)),
+        "u_turn": (Port(50, 100, 0), Port(30, 100, 180)),
+        "loop": (Port(150, 50, 90), Port(150, 40, 90)),
+        "zig_zag": (Port(20, 20, 0), Port(180, 180, 0)),
     }
-    net_widths = dict.fromkeys(netlist, 2.0)
+    net_widths = {nid: 2.0 for nid in netlist}
 
     # 3. Route
+    print("Routing complex orientation nets...")
     results = pf.route_all(netlist, net_widths)
 
     # 4. Check Results
     for nid, res in results.items():
         status = "Success" if res.is_valid else "Failed"
-        total_len = sum(comp.length for comp in res.path) if res.path else 0
-        print(f"  {nid:12}: {status}, total_length={total_len:.1f}")
+        print(f"  {nid}: {status}")
 
     # 5. Visualize
     fig, ax = plot_routing_results(results, [], bounds, netlist=netlist)

examples/06_bend_collision_models.py

@@ -30,29 +30,29 @@ def main() -> None:
     danger_map.precompute(obstacles)
 
     # We'll run three separate routers since collision_type is a router-level config
-    evaluator = CostEvaluator(engine, danger_map)
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
 
     # Scenario 1: Standard 'arc' model (High fidelity)
-    router_arc = AStarRouter(evaluator, bend_collision_type="arc")
+    router_arc = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0], bend_collision_type="arc")
     netlist_arc = {"arc_model": (Port(10, 120, 0), Port(90, 140, 90))}
 
     # Scenario 2: 'bbox' model (Conservative axis-aligned box)
-    router_bbox = AStarRouter(evaluator, bend_collision_type="bbox")
+    router_bbox = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0], bend_collision_type="bbox")
     netlist_bbox = {"bbox_model": (Port(10, 70, 0), Port(90, 90, 90))}
 
     # Scenario 3: 'clipped_bbox' model (Balanced)
-    router_clipped = AStarRouter(evaluator, bend_collision_type="clipped_bbox", bend_clip_margin=1.0)
+    router_clipped = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0], bend_collision_type="clipped_bbox", bend_clip_margin=1.0)
     netlist_clipped = {"clipped_model": (Port(10, 20, 0), Port(90, 40, 90))}
 
     # 2. Route each scenario
     print("Routing Scenario 1 (Arc)...")
-    res_arc = PathFinder(router_arc, evaluator).route_all(netlist_arc, {"arc_model": 2.0})
+    res_arc = PathFinder(router_arc, evaluator, use_tiered_strategy=False).route_all(netlist_arc, {"arc_model": 2.0})
 
     print("Routing Scenario 2 (BBox)...")
-    res_bbox = PathFinder(router_bbox, evaluator).route_all(netlist_bbox, {"bbox_model": 2.0})
+    res_bbox = PathFinder(router_bbox, evaluator, use_tiered_strategy=False).route_all(netlist_bbox, {"bbox_model": 2.0})
 
     print("Routing Scenario 3 (Clipped BBox)...")
-    res_clipped = PathFinder(router_clipped, evaluator).route_all(netlist_clipped, {"clipped_model": 2.0})
+    res_clipped = PathFinder(router_clipped, evaluator, use_tiered_strategy=False).route_all(netlist_clipped, {"clipped_model": 2.0})
 
     # 3. Combine results for visualization
     all_results = {**res_arc, **res_bbox, **res_clipped}

examples/07_large_scale_routing.py

@@ -0,0 +1,183 @@
+import numpy as np
+import time
+from inire.geometry.collision import CollisionEngine
+from inire.geometry.primitives import Port
+from inire.router.astar import AStarRouter
+from inire.router.cost import CostEvaluator
+from inire.router.danger_map import DangerMap
+from inire.router.pathfinder import PathFinder
+from inire.utils.visualization import plot_routing_results, plot_danger_map, plot_expanded_nodes, plot_expansion_density
+from shapely.geometry import box
+
+def main() -> None:
+    print("Running Example 07: Fan-Out (10 Nets, 50um Radius, 5um Grid)...")
+
+    # 1. Setup Environment
+    bounds = (0, 0, 1000, 1000)
+    engine = CollisionEngine(clearance=6.0)
+
+    # Bottleneck at x=500, 200um gap
+    obstacles = [
+        box(450, 0, 550, 400),
+        box(450, 600, 550, 1000),
+    ]
+    for obs in obstacles:
+        engine.add_static_obstacle(obs)
+
+    danger_map = DangerMap(bounds=bounds)
+    danger_map.precompute(obstacles)
+
+    evaluator = CostEvaluator(engine, danger_map, greedy_h_weight=2.0, unit_length_cost=0.1, bend_penalty=100.0, sbend_penalty=400.0, congestion_penalty=20.0)
+
+    router = AStarRouter(evaluator, node_limit=2000000, snap_size=5.0, bend_radii=[50.0], sbend_radii=[50.0], use_analytical_sbends=False)
+    pf = PathFinder(router, evaluator, max_iterations=15, base_congestion_penalty=20.0, congestion_multiplier=1.2)
+
+    # 2. Define Netlist
+    netlist = {}
+    num_nets = 10
+    start_x = 50
+    start_y_base = 500 - (num_nets * 10.0) / 2.0
+
+    end_x = 950
+    end_y_base = 100
+    end_y_pitch = 800.0 / (num_nets - 1)
+
+    for i in range(num_nets):
+        sy = round((start_y_base + i * 10.0) / 5.0) * 5.0
+        ey = round((end_y_base + i * end_y_pitch) / 5.0) * 5.0
+        netlist[f"net_{i:02d}"] = (Port(start_x, sy, 0), Port(end_x, ey, 0))
+
+    net_widths = {nid: 2.0 for nid in netlist}
+
+    # 3. Route
+    print(f"Routing {len(netlist)} nets through 200um bottleneck...")
+
+    iteration_stats = []
+
+    def iteration_callback(idx, current_results):
+        successes = sum(1 for r in current_results.values() if r.is_valid)
+        total_collisions = sum(r.collisions for r in current_results.values())
+        total_nodes = pf.router.metrics['nodes_expanded']
+
+        # Identify Hotspots
+        hotspots = {}
+        overlap_matrix = {} # (net_a, net_b) -> count
+
+        for nid, res in current_results.items():
+            if res.path:
+                for comp in res.path:
+                    for poly in comp.geometry:
+                        # Check what it overlaps with
+                        overlaps = engine.dynamic_index.intersection(poly.bounds)
+                        for other_obj_id in overlaps:
+                            other_nid, other_poly = engine.dynamic_geometries[other_obj_id]
+                            if other_nid != nid:
+                                if poly.intersects(other_poly):
+                                    # Record hotspot
+                                    cx, cy = poly.centroid.x, poly.centroid.y
+                                    grid_key = (int(cx/20)*20, int(cy/20)*20)
+                                    hotspots[grid_key] = hotspots.get(grid_key, 0) + 1
+
+                                    # Record pair
+                                    pair = tuple(sorted((nid, other_nid)))
+                                    overlap_matrix[pair] = overlap_matrix.get(pair, 0) + 1
+
+        print(f"  Iteration {idx} finished. Successes: {successes}/{len(netlist)}, Collisions: {total_collisions}")
+        if overlap_matrix:
+            top_pairs = sorted(overlap_matrix.items(), key=lambda x: x[1], reverse=True)[:3]
+            print(f"  Top Conflicts: {top_pairs}")
+        if hotspots:
+            top_hotspots = sorted(hotspots.items(), key=lambda x: x[1], reverse=True)[:3]
+            print(f"  Top Hotspots: {top_hotspots}")
+
+        # Adaptive Greediness: Decay from 2.0 to 1.1 over 25 iterations
+        new_greedy = max(1.1, 2.0 - ((idx + 1) / 25.0))
+        evaluator.greedy_h_weight = new_greedy
+        print(f"  Adaptive Greedy Weight for Next Iteration: {new_greedy:.3f}")
+
+        iteration_stats.append({
+            'Iteration': idx,
+            'Success': successes,
+            'Congestion': total_collisions,
+            'Nodes': total_nodes
+        })
+
+        # Save plots only for certain iterations to save time
+        #if idx % 20 == 0 or idx == pf.max_iterations - 1:
+        if True:
+            # Save a plot of this iteration's result
+            fig, ax = plot_routing_results(current_results, obstacles, bounds, netlist=netlist)
+            plot_danger_map(danger_map, ax=ax)
+
+            # Overlay failures: show where they stopped
+            for nid, res in current_results.items():
+                if not res.is_valid and res.path:
+                    last_p = res.path[-1].end_port
+                    target_p = netlist[nid][1]
+                    dist = abs(last_p.x - target_p.x) + abs(last_p.y - target_p.y)
+                    ax.scatter(last_p.x, last_p.y, color='red', marker='x', s=100)
+                    ax.text(last_p.x, last_p.y, f" {nid} (rem: {dist:.0f}um)", color='red', fontsize=8)
+
+            fig.savefig(f"examples/07_iteration_{idx:02d}.png")
+            import matplotlib.pyplot as plt
+            plt.close(fig)
+
+            # Plot Expansion Density if data is available
+            if pf.accumulated_expanded_nodes:
+                 fig_d, ax_d = plot_expansion_density(pf.accumulated_expanded_nodes, bounds)
+                 fig_d.savefig(f"examples/07_iteration_{idx:02d}_density.png")
+                 plt.close(fig_d)
+
+        pf.router.reset_metrics()
+
+    import cProfile, pstats
+    profiler = cProfile.Profile()
+    profiler.enable()
+    t0 = time.perf_counter()
+    results = pf.route_all(netlist, net_widths, store_expanded=True, iteration_callback=iteration_callback, shuffle_nets=True, seed=42)
+    t1 = time.perf_counter()
+    profiler.disable()
+
+    # ... (rest of the code)
+    stats = pstats.Stats(profiler).sort_stats('tottime')
+    stats.print_stats(20)
+    print(f"Routing took {t1-t0:.4f}s")
+
+    # 4. Check Results
+    print("\n--- Iteration Summary ---")
+    print(f"{'Iter':<5} | {'Success':<8} | {'Congest':<8} | {'Nodes':<10}")
+    print("-" * 40)
+    for s in iteration_stats:
+        print(f"{s['Iteration']:<5} | {s['Success']:<8} | {s['Congestion']:<8} | {s['Nodes']:<10}")
+
+    success_count = sum(1 for res in results.values() if res.is_valid)
+    print(f"\nFinal: Routed {success_count}/{len(netlist)} nets successfully.")
+
+    for nid, res in results.items():
+        target_p = netlist[nid][1]
+        if not res.is_valid:
+            last_p = res.path[-1].end_port if res.path else netlist[nid][0]
+            dist = abs(last_p.x - target_p.x) + abs(last_p.y - target_p.y)
+            print(f"  FAILED: {nid} (Stopped {dist:.1f}um from target)")
+        else:
+            types = [move.move_type for move in res.path]
+            from collections import Counter
+            counts = Counter(types)
+            print(f"  {nid}: {len(res.path)} segments, {dict(counts)}")
+
+    # 5. Visualize
+    fig, ax = plot_routing_results(results, obstacles, bounds, netlist=netlist)
+
+    # Overlay Danger Map
+    plot_danger_map(danger_map, ax=ax)
+
+    # Overlay Expanded Nodes from last routed net (as an example)
+    if pf.router.last_expanded_nodes:
+        print(f"Plotting {len(pf.router.last_expanded_nodes)} expanded nodes for the last net...")
+        plot_expanded_nodes(pf.router.last_expanded_nodes, ax=ax, color='blue', alpha=0.1)
+
+    fig.savefig("examples/07_large_scale_routing.png")
+    print("Saved plot to examples/07_large_scale_routing.png")
+
+if __name__ == "__main__":
+    main()

examples/08_custom_bend_geometry.py

@@ -0,0 +1,54 @@
+from shapely.geometry import Polygon
+
+from inire.geometry.collision import CollisionEngine
+from inire.geometry.primitives import Port
+from inire.router.astar import AStarRouter
+from inire.router.cost import CostEvaluator
+from inire.router.danger_map import DangerMap
+from inire.router.pathfinder import PathFinder
+from inire.utils.visualization import plot_routing_results
+
+
+def main() -> None:
+    print("Running Example 08: Custom Bend Geometry...")
+
+    # 1. Setup Environment
+    bounds = (0, 0, 150, 150)
+    engine = CollisionEngine(clearance=2.0)
+    danger_map = DangerMap(bounds=bounds)
+    danger_map.precompute([])
+
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    router = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
+    pf = PathFinder(router, evaluator)
+
+    # 2. Define Netlist
+    netlist = {
+        "custom_bend": (Port(20, 20, 0), Port(100, 100, 90)),
+    }
+    net_widths = {"custom_bend": 2.0}
+
+    # 3. Route with standard arc first
+    print("Routing with standard arc...")
+    results_std = pf.route_all(netlist, net_widths)
+
+    # 4. Define a custom 'trapezoid' bend model
+    # (Just for demonstration - we override the collision model during search)
+    custom_poly = Polygon([(0, 0), (20, 0), (20, 20), (0, 20)]) # Oversized box
+    
+    print("Routing with custom collision model...")
+    # Override bend_collision_type with a literal Polygon
+    router_custom = AStarRouter(evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0], bend_collision_type=custom_poly)
+    results_custom = PathFinder(router_custom, evaluator, use_tiered_strategy=False).route_all(
+        {"custom_model": netlist["custom_bend"]}, {"custom_model": 2.0}
+    )
+
+    # 5. Visualize
+    all_results = {**results_std, **results_custom}
+    fig, ax = plot_routing_results(all_results, [], bounds, netlist=netlist)
+    fig.savefig("examples/08_custom_bend_geometry.png")
+    print("Saved plot to examples/08_custom_bend_geometry.png")
+
+
+if __name__ == "__main__":
+    main()

examples/09_unroutable_best_effort.py

@@ -0,0 +1,59 @@
+from inire.geometry.collision import CollisionEngine
+from inire.geometry.primitives import Port
+from inire.router.astar import AStarRouter
+from inire.router.cost import CostEvaluator
+from inire.router.danger_map import DangerMap
+from inire.router.pathfinder import PathFinder
+from inire.utils.visualization import plot_routing_results
+from shapely.geometry import box
+
+def main() -> None:
+    print("Running Example 09: Best-Effort (Unroutable Net)...")
+
+    # 1. Setup Environment
+    bounds = (0, 0, 100, 100)
+    engine = CollisionEngine(clearance=2.0)
+    
+    # Create a 'cage' that completely blocks the target
+    cage = [
+        box(70, 30, 75, 70), # Left wall
+        box(70, 70, 95, 75), # Top wall
+        box(70, 25, 95, 30), # Bottom wall
+    ]
+    for obs in cage:
+        engine.add_static_obstacle(obs)
+
+    danger_map = DangerMap(bounds=bounds)
+    danger_map.precompute(cage)
+
+    evaluator = CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
+    # Use a low node limit to fail faster
+    router = AStarRouter(evaluator, node_limit=2000, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
+    
+    # Enable partial path return
+    pf = PathFinder(router, evaluator)
+
+    # 2. Define Netlist: start outside, target inside the cage
+    netlist = {
+        "trapped_net": (Port(10, 50, 0), Port(85, 50, 0)),
+    }
+    net_widths = {"trapped_net": 2.0}
+
+    # 3. Route
+    print("Routing net into a cage (should fail and return partial)...")
+    results = pf.route_all(netlist, net_widths)
+
+    # 4. Check Results
+    res = results["trapped_net"]
+    if not res.is_valid:
+        print(f"Net failed to route as expected. Partial path length: {len(res.path)} segments.")
+    else:
+        print("Wait, it found a way in? Check the cage geometry!")
+
+    # 5. Visualize
+    fig, ax = plot_routing_results(results, cage, bounds, netlist=netlist)
+    fig.savefig("examples/09_unroutable_best_effort.png")
+    print("Saved plot to examples/09_unroutable_best_effort.png")
+
+if __name__ == "__main__":
+    main()

examples/README.md

@@ -0,0 +1,38 @@
+# Inire Routing Examples
+
+This directory contains examples demonstrating the features and architectural capabilities of the `inire` router.
+
+## Architectural Visualization
+In all plots generated by `inire`, we distinguish between the search-time geometry and the final "actual" geometry:
+*   **Dashed Lines & Translucent Fill**: The **Collision Proxy** used during the A* search (e.g., `clipped_bbox` or `bbox`). This represents the conservative envelope the router used to guarantee clearance.
+*   **Solid Lines**: The **Actual Geometry** (high-fidelity arcs). This is the exact shape that will be used for PDK generation and fabrication.
+
+---
+
+## 1. Fan-Out (Negotiated Congestion)
+Demonstrates the Negotiated Congestion algorithm handling multiple intersecting nets. The router iteratively increases penalties for overlaps until a collision-free solution is found. This example shows a bundle of nets fanning out through a narrow bottleneck.
+
+![Fan-Out Routing](07_large_scale_routing.png)
+
+## 2. Custom Bend Geometry Models
+`inire` supports multiple collision models for bends, allowing a trade-off between search speed and geometric accuracy:
+*   **Arc**: High-fidelity geometry (Highest accuracy).
+*   **BBox**: Simple axis-aligned bounding box (Fastest search).
+*   **Clipped BBox**: A balanced model that clips the corners of the AABB to better fit the arc (Optimal performance).
+
+![Custom Bend Geometry](08_custom_bend_geometry.png)
+
+## 3. Unroutable Nets & Best-Effort Display
+When a net is physically blocked or exceeds the node limit, the router returns the "best-effort" partial path—the path that reached the point closest to the target according to the heuristic. This is critical for debugging design constraints.
+
+![Best Effort Display](09_unroutable_best_effort.png)
+
+## 4. Orientation Stress Test
+Demonstrates the router's ability to handle complex orientation requirements, including U-turns, 90-degree flips, and loops.
+
+![Orientation Stress Test](05_orientation_stress.png)
+
+## 5. Tiered Fidelity & Lazy Dilation
+Our architecture leverages two key optimizations for high-performance routing:
+1.  **Tiered Fidelity**: Initial routing passes use fast `clipped_bbox` proxies. If collisions are found, the system automatically escalates to high-fidelity `arc` geometry for the affected regions.
+2.  **Lazy Dilation**: Geometric buffering (dilation) is deferred until a collision check is strictly necessary, avoiding thousands of redundant `buffer()` and `translate()` calls.

inire/geometry/collision.py

@@ -2,12 +2,16 @@ from __future__ import annotations
 
 from typing import TYPE_CHECKING, Literal
 import rtree
+import numpy
 from shapely.prepared import prep
+from shapely.strtree import STRtree
+from shapely.geometry import box
 
 if TYPE_CHECKING:
     from shapely.geometry import Polygon
     from shapely.prepared import PreparedGeometry
     from inire.geometry.primitives import Port
+    from inire.geometry.components import ComponentResult
 
 
 class CollisionEngine:
@@ -16,8 +20,12 @@ class CollisionEngine:
     """
     __slots__ = (
         'clearance', 'max_net_width', 'safety_zone_radius',
-        'static_index', 'static_geometries', 'static_dilated', 'static_prepared', '_static_id_counter',
-        'dynamic_index', 'dynamic_geometries', 'dynamic_dilated', 'dynamic_prepared', '_dynamic_id_counter'
+        'static_index', 'static_geometries', 'static_dilated', 'static_prepared', 
+        'static_is_rect', 'static_tree', 'static_obj_ids', 'static_safe_cache', 
+        'static_grid', 'grid_cell_size', '_static_id_counter',
+        'dynamic_index', 'dynamic_geometries', 'dynamic_dilated', 'dynamic_prepared', 
+        'dynamic_tree', 'dynamic_obj_ids', 'dynamic_grid', '_dynamic_id_counter',
+        'metrics'
     )
 
     clearance: float
@@ -52,16 +60,54 @@ class CollisionEngine:
         self.static_geometries: dict[int, Polygon] = {}  # ID -> Raw Polygon
         self.static_dilated: dict[int, Polygon] = {}     # ID -> Dilated Polygon (by clearance)
         self.static_prepared: dict[int, PreparedGeometry] = {}  # ID -> Prepared Dilated
+        self.static_is_rect: dict[int, bool] = {} # Optimization for ray_cast
+        self.static_tree: STRtree | None = None
+        self.static_obj_ids: list[int] = [] # Mapping from tree index to obj_id
+        self.static_safe_cache: set[tuple] = set() # Global cache for safe move-port combinations
+        self.static_grid: dict[tuple[int, int], list[int]] = {}
+        self.grid_cell_size = 50.0 # 50um grid cells for broad phase
         self._static_id_counter = 0
 
         # Dynamic paths for multi-net congestion
         self.dynamic_index = rtree.index.Index()
-        # obj_id -> (net_id, raw_geometry)
         self.dynamic_geometries: dict[int, tuple[str, Polygon]] = {}
-        # obj_id -> dilated_geometry (by clearance/2)
         self.dynamic_dilated: dict[int, Polygon] = {}
         self.dynamic_prepared: dict[int, PreparedGeometry] = {}
+        self.dynamic_tree: STRtree | None = None
+        self.dynamic_obj_ids: list[int] = []
+        self.dynamic_grid: dict[tuple[int, int], list[int]] = {}
         self._dynamic_id_counter = 0
+        
+        self.metrics = {
+            'static_cache_hits': 0,
+            'static_grid_skips': 0,
+            'static_tree_queries': 0,
+            'static_straight_fast': 0,
+            'congestion_grid_skips': 0,
+            'congestion_tree_queries': 0,
+            'safety_zone_checks': 0
+        }
+
+    def reset_metrics(self) -> None:
+        """ Reset all performance counters. """
+        for k in self.metrics:
+            self.metrics[k] = 0
+
+    def get_metrics_summary(self) -> str:
+        """ Return a human-readable summary of collision performance. """
+        m = self.metrics
+        total_static = m['static_cache_hits'] + m['static_grid_skips'] + m['static_tree_queries'] + m['static_straight_fast']
+        static_eff = ((m['static_cache_hits'] + m['static_grid_skips'] + m['static_straight_fast']) / total_static * 100) if total_static > 0 else 0
+        
+        total_cong = m['congestion_grid_skips'] + m['congestion_tree_queries']
+        cong_eff = (m['congestion_grid_skips'] / total_cong * 100) if total_cong > 0 else 0
+        
+        return (f"Collision Performance: \n"
+                f"  Static: {total_static} checks, {static_eff:.1f}% bypassed STRtree\n"
+                f"    (Cache={m['static_cache_hits']}, Grid={m['static_grid_skips']}, StraightFast={m['static_straight_fast']}, Tree={m['static_tree_queries']})\n"
+                f"  Congestion: {total_cong} checks, {cong_eff:.1f}% bypassed STRtree\n"
+                f"    (Grid={m['congestion_grid_skips']}, Tree={m['congestion_tree_queries']})\n"
+                f"  Safety Zone: {m['safety_zone_checks']} full intersections performed")
 
     def add_static_obstacle(self, polygon: Polygon) -> None:
         """
@@ -73,11 +119,46 @@ class CollisionEngine:
         obj_id = self._static_id_counter
         self._static_id_counter += 1
 
-        dilated = polygon.buffer(self.clearance)
+        # Use MITRE join style to preserve rectangularity of boxes
+        dilated = polygon.buffer(self.clearance, join_style=2)
         self.static_geometries[obj_id] = polygon
         self.static_dilated[obj_id] = dilated
         self.static_prepared[obj_id] = prep(dilated)
         self.static_index.insert(obj_id, dilated.bounds)
+        
+        # Invalidate higher-level spatial data
+        self.static_tree = None
+        self.static_grid = {} # Rebuild on demand
+        
+        # Check if it's an axis-aligned rectangle (approximately)
+        # Dilated rectangle of an axis-aligned rectangle IS an axis-aligned rectangle.
+        b = dilated.bounds
+        area = (b[2] - b[0]) * (b[3] - b[1])
+        if abs(dilated.area - area) < 1e-4:
+             self.static_is_rect[obj_id] = True
+        else:
+             self.static_is_rect[obj_id] = False
+
+    def _ensure_static_tree(self) -> None:
+        if self.static_tree is None and self.static_dilated:
+            ids = sorted(self.static_dilated.keys())
+            geoms = [self.static_dilated[i] for i in ids]
+            self.static_tree = STRtree(geoms)
+            self.static_obj_ids = ids
+
+    def _ensure_static_grid(self) -> None:
+        if not self.static_grid and self.static_dilated:
+            cs = self.grid_cell_size
+            for obj_id, poly in self.static_dilated.items():
+                b = poly.bounds
+                min_gx, max_gx = int(b[0] / cs), int(b[2] / cs)
+                min_gy, max_gy = int(b[1] / cs), int(b[3] / cs)
+                for gx in range(min_gx, max_gx + 1):
+                    for gy in range(min_gy, max_gy + 1):
+                        cell = (gx, gy)
+                        if cell not in self.static_grid:
+                            self.static_grid[cell] = []
+                        self.static_grid[cell].append(obj_id)
 
     def add_path(self, net_id: str, geometry: list[Polygon], dilated_geometry: list[Polygon] | None = None) -> None:
         """
@@ -99,6 +180,30 @@ class CollisionEngine:
             self.dynamic_dilated[obj_id] = dil
             self.dynamic_prepared[obj_id] = prep(dil)
             self.dynamic_index.insert(obj_id, dil.bounds)
+        
+        self.dynamic_tree = None
+        self.dynamic_grid = {}
+
+    def _ensure_dynamic_tree(self) -> None:
+        if self.dynamic_tree is None and self.dynamic_dilated:
+            ids = sorted(self.dynamic_dilated.keys())
+            geoms = [self.dynamic_dilated[i] for i in ids]
+            self.dynamic_tree = STRtree(geoms)
+            self.dynamic_obj_ids = ids
+
+    def _ensure_dynamic_grid(self) -> None:
+        if not self.dynamic_grid and self.dynamic_dilated:
+            cs = self.grid_cell_size
+            for obj_id, poly in self.dynamic_dilated.items():
+                b = poly.bounds
+                min_gx, max_gx = int(b[0] / cs), int(b[2] / cs)
+                min_gy, max_gy = int(b[1] / cs), int(b[3] / cs)
+                for gx in range(min_gx, max_gx + 1):
+                    for gy in range(min_gy, max_gy + 1):
+                        cell = (gx, gy)
+                        if cell not in self.dynamic_grid:
+                            self.dynamic_grid[cell] = []
+                        self.dynamic_grid[cell].append(obj_id)
 
     def remove_path(self, net_id: str) -> None:
         """
@@ -113,6 +218,10 @@ class CollisionEngine:
             dilated = self.dynamic_dilated.pop(obj_id)
             self.dynamic_prepared.pop(obj_id)
             self.dynamic_index.delete(obj_id, dilated.bounds)
+        
+        if to_remove:
+            self.dynamic_tree = None
+            self.dynamic_grid = {}
 
     def lock_net(self, net_id: str) -> None:
         """
@@ -152,6 +261,279 @@ class CollisionEngine:
         res = self.check_collision(geometry, net_id, buffer_mode='congestion')
         return int(res)
 
+    def check_move_straight_static(
+            self,
+            origin: Port,
+            length: float,
+            ) -> bool:
+        """
+        Specialized fast static check for Straights.
+        """
+        self.metrics['static_straight_fast'] += 1
+        # FAST PATH: Grid check
+        self._ensure_static_grid()
+        cs = self.grid_cell_size
+        
+        rad = numpy.radians(origin.orientation)
+        dx = length * numpy.cos(rad)
+        dy = length * numpy.sin(rad)
+        
+        # Move bounds
+        xmin, xmax = sorted([origin.x, origin.x + dx])
+        ymin, ymax = sorted([origin.y, origin.y + dy])
+        
+        # Inflate by clearance/2 for waveguide half-width? 
+        # No, static obstacles are ALREADY inflated by full clearance.
+        # So we just check if the centerline hits an inflated obstacle.
+        
+        min_gx, max_gx = int(xmin / cs), int(xmax / cs)
+        min_gy, max_gy = int(ymin / cs), int(ymax / cs)
+        
+        static_grid = self.static_grid
+        static_dilated = self.static_dilated
+        static_is_rect = self.static_is_rect
+        static_prepared = self.static_prepared
+        
+        inv_dx = 1.0/dx if abs(dx) > 1e-12 else 1e30
+        inv_dy = 1.0/dy if abs(dy) > 1e-12 else 1e30
+
+        checked_ids = set()
+        for gx in range(min_gx, max_gx + 1):
+            for gy in range(min_gy, max_gy + 1):
+                if (gx, gy) in static_grid:
+                    for obj_id in static_grid[(gx, gy)]:
+                        if obj_id in checked_ids: continue
+                        checked_ids.add(obj_id)
+                        
+                        b = static_dilated[obj_id].bounds
+                        # Slab Method
+                        if abs(dx) < 1e-12:
+                            if origin.x < b[0] or origin.x > b[2]: continue
+                            tx_min, tx_max = -1e30, 1e30
+                        else:
+                            tx_min = (b[0] - origin.x) * inv_dx
+                            tx_max = (b[2] - origin.x) * inv_dx
+                            if tx_min > tx_max: tx_min, tx_max = tx_max, tx_min
+                            
+                        if abs(dy) < 1e-12:
+                            if origin.y < b[1] or origin.y > b[3]: continue
+                            ty_min, ty_max = -1e30, 1e30
+                        else:
+                            ty_min = (b[1] - origin.y) * inv_dy
+                            ty_max = (b[3] - origin.y) * inv_dy
+                            if ty_min > ty_max: ty_min, ty_max = ty_max, ty_min
+                            
+                        t_min = max(tx_min, ty_min)
+                        t_max = min(tx_max, ty_max)
+                        
+                        if t_max <= 1e-9 or t_min > t_max or t_min >= 1.0 - 1e-9:
+                            continue
+                            
+                        # If rectangle, slab is exact
+                        if static_is_rect[obj_id]:
+                            return True
+                            
+                        # Fallback for complex obstacles
+                        # (We could still use ray_cast here but we want exact)
+                        # For now, if hits AABB, check prepared
+                        from shapely.geometry import LineString
+                        line = LineString([(origin.x, origin.y), (origin.x+dx, origin.y+dy)])
+                        if static_prepared[obj_id].intersects(line):
+                            return True
+        return False
+
+    def check_move_static(
+            self,
+            result: ComponentResult,
+            start_port: Port | None = None,
+            end_port: Port | None = None,
+            ) -> bool:
+        """
+        Check if a move (ComponentResult) hits any static obstacles.
+        """
+        # FAST PATH 1: Safety cache check
+        cache_key = (result.move_type, 
+                     round(start_port.x, 3) if start_port else 0, 
+                     round(start_port.y, 3) if start_port else 0,
+                     round(end_port.x, 3) if end_port else 0, 
+                     round(end_port.y, 3) if end_port else 0)
+        
+        if cache_key in self.static_safe_cache:
+            self.metrics['static_cache_hits'] += 1
+            return False
+
+        # FAST PATH 2: Spatial grid check (bypasses STRtree for empty areas)
+        self._ensure_static_grid()
+        cs = self.grid_cell_size
+        b = result.total_bounds
+        min_gx, max_gx = int(b[0] / cs), int(b[2] / cs)
+        min_gy, max_gy = int(b[1] / cs), int(b[3] / cs)
+        
+        any_candidates = False
+        static_grid = self.static_grid
+        for gx in range(min_gx, max_gx + 1):
+            for gy in range(min_gy, max_gy + 1):
+                if (gx, gy) in static_grid:
+                    any_candidates = True
+                    break
+            if any_candidates: break
+            
+        if not any_candidates:
+            self.metrics['static_grid_skips'] += 1
+            self.static_safe_cache.add(cache_key)
+            return False
+
+        self.metrics['static_tree_queries'] += 1
+        self._ensure_static_tree()
+        if self.static_tree is None:
+            return False
+
+        # Vectorized Broad phase + Narrow phase
+        # Pass all polygons in the move at once
+        res_indices, tree_indices = self.static_tree.query(result.geometry, predicate='intersects')
+        if tree_indices.size == 0:
+            self.static_safe_cache.add(cache_key)
+            return False
+
+        # If we have hits, we must check safety zones
+        static_obj_ids = self.static_obj_ids
+        for i in range(tree_indices.size):
+            poly_idx = res_indices[i]
+            hit_idx = tree_indices[i]
+            obj_id = static_obj_ids[hit_idx]
+            poly = result.geometry[poly_idx]
+            if self._is_in_safety_zone(poly, obj_id, start_port, end_port):
+                continue
+            return True
+        
+        self.static_safe_cache.add(cache_key)
+        return False
+
+    def check_move_congestion(
+            self,
+            result: ComponentResult,
+            net_id: str,
+            ) -> int:
+        """
+        Count overlaps of a move with other dynamic paths.
+        """
+        if result.total_dilated_bounds_box is None:
+             return 0 
+
+        # FAST PATH: Grid check
+        self._ensure_dynamic_grid()
+        if not self.dynamic_grid:
+            return 0
+            
+        cs = self.grid_cell_size
+        b = result.total_dilated_bounds
+        min_gx, max_gx = int(b[0] / cs), int(b[2] / cs)
+        min_gy, max_gy = int(b[1] / cs), int(b[3] / cs)
+        
+        any_candidates = False
+        dynamic_grid = self.dynamic_grid
+        dynamic_geometries = self.dynamic_geometries
+        for gx in range(min_gx, max_gx + 1):
+            for gy in range(min_gy, max_gy + 1):
+                cell = (gx, gy)
+                if cell in dynamic_grid:
+                    # Check if any obj_id in this cell belongs to another net
+                    for obj_id in dynamic_grid[cell]:
+                        other_net_id, _ = dynamic_geometries[obj_id]
+                        if other_net_id != net_id:
+                            any_candidates = True
+                            break
+                if any_candidates: break
+            if any_candidates: break
+            
+        if not any_candidates:
+            self.metrics['congestion_grid_skips'] += 1
+            return 0
+
+        # SLOW PATH: STRtree
+        self.metrics['congestion_tree_queries'] += 1
+        self._ensure_dynamic_tree()
+        if self.dynamic_tree is None:
+            return 0
+
+        # Vectorized query: pass the whole list of polygons
+        # result.dilated_geometry is list[Polygon]
+        # query() returns (2, M) array of [geometry_indices, tree_indices]
+        res_indices, tree_indices = self.dynamic_tree.query(result.dilated_geometry, predicate='intersects')
+        if tree_indices.size == 0:
+            return 0
+
+        count = 0
+        dynamic_geometries = self.dynamic_geometries
+        dynamic_obj_ids = self.dynamic_obj_ids
+        
+        # We need to filter by net_id and count UNIQUE overlaps?
+        # Actually, if a single move polygon hits multiple other net polygons, it's multiple overlaps.
+        # But if multiple move polygons hit the SAME other net polygon, is it multiple overlaps?
+        # Usually, yes, because cost is proportional to volume of overlap.
+        for hit_idx in tree_indices:
+            obj_id = dynamic_obj_ids[hit_idx]
+            other_net_id, _ = dynamic_geometries[obj_id]
+            if other_net_id != net_id:
+                count += 1
+        return count
+
+    def _is_in_safety_zone(self, geometry: Polygon, obj_id: int, start_port: Port | None, end_port: Port | None) -> bool:
+        """ Helper to check if an intersection is within a port safety zone. """
+        sz = self.safety_zone_radius
+        static_dilated = self.static_dilated
+        
+        # Optimization: Skip expensive intersection if neither port is near the obstacle's bounds
+        is_near_port = False
+        b = static_dilated[obj_id].bounds
+        if start_port:
+            if (b[0] - sz <= start_port.x <= b[2] + sz and 
+                b[1] - sz <= start_port.y <= b[3] + sz):
+                is_near_port = True
+        if not is_near_port and end_port:
+            if (b[0] - sz <= end_port.x <= b[2] + sz and 
+                b[1] - sz <= end_port.y <= b[3] + sz):
+                is_near_port = True
+        
+        if not is_near_port:
+            return False # Collision is NOT in safety zone
+
+        # Only if near port, do the expensive check
+        self.metrics['safety_zone_checks'] += 1
+        raw_obstacle = self.static_geometries[obj_id]
+        intersection = geometry.intersection(raw_obstacle)
+        if intersection.is_empty:
+            return True # Not actually hitting the RAW obstacle (only the buffer)
+            
+        ix_bounds = intersection.bounds
+        # Check start port
+        if start_port:
+            if (abs(ix_bounds[0] - start_port.x) < sz and
+                abs(ix_bounds[2] - start_port.x) < sz and
+                abs(ix_bounds[1] - start_port.y) < sz and
+                abs(ix_bounds[3] - start_port.y) < sz):
+                return True # Is safe
+        # Check end port
+        if end_port:
+            if (abs(ix_bounds[0] - end_port.x) < sz and
+                abs(ix_bounds[2] - end_port.x) < sz and
+                abs(ix_bounds[1] - end_port.y) < sz and
+                abs(ix_bounds[3] - end_port.y) < sz):
+                return True # Is safe
+                
+        return False
+
+    def check_congestion(
+            self,
+            geometry: Polygon,
+            net_id: str,
+            dilated_geometry: Polygon | None = None,
+            ) -> int:
+        """
+        Alias for check_collision(buffer_mode='congestion') for backward compatibility.
+        """
+        res = self.check_collision(geometry, net_id, buffer_mode='congestion', dilated_geometry=dilated_geometry)
+        return int(res)
     def check_collision(
             self,
             geometry: Polygon,
@@ -160,71 +542,165 @@ class CollisionEngine:
             start_port: Port | None = None,
             end_port: Port | None = None,
             dilated_geometry: Polygon | None = None,
+            bounds: tuple[float, float, float, float] | None = None,
             ) -> bool | int:
         """
         Check for collisions using unified dilation logic.
-
-        Args:
-            geometry: Raw geometry to check.
-            net_id: Identifier for the net.
-            buffer_mode: 'static' (full clearance) or 'congestion' (shared).
-            start_port: Optional start port for safety zone.
-            end_port: Optional end port for safety zone.
-            dilated_geometry: Optional pre-buffered geometry (clearance/2).
-
-        Returns:
-            Boolean if static, integer count if congestion.
         """
         if buffer_mode == 'static':
-            # Use raw query against pre-dilated obstacles
-            candidates = self.static_index.intersection(geometry.bounds)
-
-            for obj_id in candidates:
-                if self.static_prepared[obj_id].intersects(geometry):
-                    if start_port or end_port:
-                        # Optimization: Skip expensive intersection if neither port is near the obstacle's bounds
-                        # (Plus a small margin for safety zone)
-                        sz = self.safety_zone_radius
-                        is_near_port = False
-                        for p in [start_port, end_port]:
-                            if p:
-                                # Quick bounds check
-                                b = self.static_dilated[obj_id].bounds
-                                if (b[0] - sz <= p.x <= b[2] + sz and 
-                                    b[1] - sz <= p.y <= b[3] + sz):
-                                    is_near_port = True
-                                    break
-                        
-                        if not is_near_port:
-                            return True # Collision, and not near any port safety zone
-
-                        # Only if near port, do the expensive check
-                        raw_obstacle = self.static_geometries[obj_id]
-                        intersection = geometry.intersection(raw_obstacle)
-                        if not intersection.is_empty:
-                            ix_minx, ix_miny, ix_maxx, ix_maxy = intersection.bounds
-                            is_safe = False
-                            for p in [start_port, end_port]:
-                                if p and (abs(ix_minx - p.x) < sz and
-                                          abs(ix_maxx - p.x) < sz and
-                                          abs(ix_miny - p.y) < sz and
-                                          abs(ix_maxy - p.y) < sz):
-                                    is_safe = True
-                                    break
-                            
-                            if is_safe:
-                                continue
-                    return True
+            self._ensure_static_tree()
+            if self.static_tree is None:
+                return False
+                
+            hits = self.static_tree.query(geometry, predicate='intersects')
+            static_obj_ids = self.static_obj_ids
+            for hit_idx in hits:
+                obj_id = static_obj_ids[hit_idx]
+                if self._is_in_safety_zone(geometry, obj_id, start_port, end_port):
+                    continue
+                return True
             return False
 
         # buffer_mode == 'congestion'
+        self._ensure_dynamic_tree()
+        if self.dynamic_tree is None:
+            return 0
+            
         dilation = self.clearance / 2.0
         test_poly = dilated_geometry if dilated_geometry else geometry.buffer(dilation)
-        candidates = self.dynamic_index.intersection(test_poly.bounds)
-
+        
+        hits = self.dynamic_tree.query(test_poly, predicate='intersects')
         count = 0
-        for obj_id in candidates:
-            other_net_id, _ = self.dynamic_geometries[obj_id]
-            if other_net_id != net_id and self.dynamic_prepared[obj_id].intersects(test_poly):
+        dynamic_geometries = self.dynamic_geometries
+        dynamic_obj_ids = self.dynamic_obj_ids
+        
+        for hit_idx in hits:
+            obj_id = dynamic_obj_ids[hit_idx]
+            other_net_id, _ = dynamic_geometries[obj_id]
+            if other_net_id != net_id:
                 count += 1
         return count
+
+    def ray_cast(self, origin: Port, angle_deg: float, max_dist: float = 2000.0) -> float:
+        """
+        Cast a ray and find the distance to the nearest static obstacle.
+
+        Args:
+            origin: Starting port (x, y).
+            angle_deg: Ray direction in degrees.
+            max_dist: Maximum lookahead distance.
+
+        Returns:
+            Distance to first collision, or max_dist if clear.
+        """
+        import numpy
+        from shapely.geometry import LineString
+
+        rad = numpy.radians(angle_deg)
+        cos_val = numpy.cos(rad)
+        sin_val = numpy.sin(rad)
+        dx = max_dist * cos_val
+        dy = max_dist * sin_val
+        
+        # 1. Pre-calculate ray direction inverses for fast slab intersection
+        # Use a small epsilon to avoid divide by zero, but handle zero dx/dy properly.
+        if abs(dx) < 1e-12:
+            inv_dx = 1e30 # Represent infinity
+        else:
+            inv_dx = 1.0 / dx
+            
+        if abs(dy) < 1e-12:
+            inv_dy = 1e30 # Represent infinity
+        else:
+            inv_dy = 1.0 / dy
+        
+        # Ray AABB for initial R-Tree query
+        min_x, max_x = sorted([origin.x, origin.x + dx])
+        min_y, max_y = sorted([origin.y, origin.y + dy])
+        
+        # 1. Query R-Tree
+        candidates = list(self.static_index.intersection((min_x, min_y, max_x, max_y)))
+        if not candidates:
+            return max_dist
+        
+        min_dist = max_dist
+        
+        # 2. Check Intersections
+        # Note: We intersect with DILATED obstacles to account for clearance
+        static_dilated = self.static_dilated
+        static_prepared = self.static_prepared
+        
+        # Optimization: Sort candidates by approximate distance to origin
+        # (Using a simpler distance measure for speed)
+        def approx_dist_sq(obj_id):
+            b = static_dilated[obj_id].bounds
+            return (b[0] - origin.x)**2 + (b[1] - origin.y)**2
+            
+        candidates.sort(key=approx_dist_sq)
+
+        ray_line = None # Lazy creation
+
+        for obj_id in candidates:
+            b = static_dilated[obj_id].bounds
+            
+            # Fast Ray-Box intersection (Slab Method)
+            # Correctly handle potential for dx=0 or dy=0
+            if abs(dx) < 1e-12:
+                if origin.x < b[0] or origin.x > b[2]:
+                    continue
+                tx_min, tx_max = -1e30, 1e30
+            else:
+                tx_min = (b[0] - origin.x) * inv_dx
+                tx_max = (b[2] - origin.x) * inv_dx
+                if tx_min > tx_max: tx_min, tx_max = tx_max, tx_min
+            
+            if abs(dy) < 1e-12:
+                if origin.y < b[1] or origin.y > b[3]:
+                    continue
+                ty_min, ty_max = -1e30, 1e30
+            else:
+                ty_min = (b[1] - origin.y) * inv_dy
+                ty_max = (b[3] - origin.y) * inv_dy
+                if ty_min > ty_max: ty_min, ty_max = ty_max, ty_min
+            
+            t_min = max(tx_min, ty_min)
+            t_max = min(tx_max, ty_max)
+            
+            # Intersection if [t_min, t_max] intersects [0, 1]
+            if t_max < 0 or t_min > t_max or t_min >= (min_dist / max_dist) or t_min > 1.0:
+                continue
+            
+            # Optimization: If it's a rectangle, the slab result is exact!
+            if self.static_is_rect[obj_id]:
+                min_dist = max(0.0, t_min * max_dist)
+                continue
+
+            # If we are here, the ray hits the AABB. Now check the actual polygon.
+            if ray_line is None:
+                ray_line = LineString([(origin.x, origin.y), (origin.x + dx, origin.y + dy)])
+
+            if static_prepared[obj_id].intersects(ray_line):
+                # Calculate exact intersection distance
+                intersection = ray_line.intersection(static_dilated[obj_id])
+                if intersection.is_empty:
+                    continue
+                
+                # Intersection could be MultiLineString or LineString or Point
+                def get_dist(geom):
+                    if hasattr(geom, 'geoms'): # Multi-part
+                        return min(get_dist(g) for g in geom.geoms)
+                    # For line string, the intersection is the segment INSIDE the obstacle.
+                    coords = geom.coords
+                    p1 = coords[0]
+                    return numpy.sqrt((p1[0] - origin.x)**2 + (p1[1] - origin.y)**2)
+
+                try:
+                    d = get_dist(intersection)
+                    if d < min_dist:
+                        min_dist = d
+                        # Update ray_line for more aggressive pruning? 
+                        # Actually just update min_dist and we use it in the t_min check.
+                except Exception:
+                    pass
+
+        return min_dist

inire/geometry/components.py

@@ -1,55 +1,36 @@
 from __future__ import annotations
 
-from typing import Literal, cast
+import math
+from typing import Literal, cast, Any
 import numpy
 import shapely
-from shapely.geometry import Polygon, box
+from shapely.geometry import Polygon, box, MultiPolygon
 from shapely.ops import unary_union
 
 from .primitives import Port
 
 
 
-# Search Grid Snap (1.0 µm)
-SEARCH_GRID_SNAP_UM = 1.0
+# Search Grid Snap (5.0 µm default)
+SEARCH_GRID_SNAP_UM = 5.0
 
 
-def snap_search_grid(value: float) -> float:
+def snap_search_grid(value: float, snap_size: float = SEARCH_GRID_SNAP_UM) -> float:
     """
     Snap a coordinate to the nearest search grid unit.
-
-    Args:
-        value: Value to snap.
-
-    Returns:
-        Snapped value.
     """
-    return round(value / SEARCH_GRID_SNAP_UM) * SEARCH_GRID_SNAP_UM
+    return round(value / snap_size) * snap_size
 
 
 class ComponentResult:
     """
-    The result of a component generation: geometry, final port, and physical length.
+    Standard container for generated move geometry and state.
     """
-    __slots__ = ('geometry', 'dilated_geometry', 'end_port', 'length', 'bounds', 'dilated_bounds')
-
-    geometry: list[Polygon]
-    """ List of polygons representing the component geometry """
-
-    dilated_geometry: list[Polygon] | None
-    """ Optional list of pre-dilated polygons for collision optimization """
-
-    end_port: Port
-    """ The final port after the component """
-
-    length: float
-    """ Physical length of the component path """
-
-    bounds: numpy.ndarray
-    """ Pre-calculated bounds for each polygon in geometry """
-
-    dilated_bounds: numpy.ndarray | None
-    """ Pre-calculated bounds for each polygon in dilated_geometry """
+    __slots__ = (
+        'geometry', 'dilated_geometry', 'proxy_geometry', 'actual_geometry',
+        'end_port', 'length', 'move_type', 'bounds', 'dilated_bounds', 
+        'total_bounds', 'total_dilated_bounds', 'total_bounds_box', 'total_dilated_bounds_box', '_t_cache'
+    )
 
     def __init__(
             self,
@@ -57,38 +38,99 @@ class ComponentResult:
             end_port: Port,
             length: float,
             dilated_geometry: list[Polygon] | None = None,
+            proxy_geometry: list[Polygon] | None = None,
+            actual_geometry: list[Polygon] | None = None,
+            skip_bounds: bool = False,
+            move_type: str = 'Unknown'
             ) -> None:
         self.geometry = geometry
         self.dilated_geometry = dilated_geometry
+        self.proxy_geometry = proxy_geometry
+        self.actual_geometry = actual_geometry
         self.end_port = end_port
         self.length = length
-        # Vectorized bounds calculation
-        self.bounds = shapely.bounds(geometry)
-        self.dilated_bounds = shapely.bounds(dilated_geometry) if dilated_geometry is not None else None
+        self.move_type = move_type
+        self._t_cache = {}
+        if not skip_bounds:
+            # Vectorized bounds calculation
+            self.bounds = shapely.bounds(geometry)
+            # Total bounds across all polygons in the move
+            self.total_bounds = numpy.array([
+                numpy.min(self.bounds[:, 0]),
+                numpy.min(self.bounds[:, 1]),
+                numpy.max(self.bounds[:, 2]),
+                numpy.max(self.bounds[:, 3])
+            ])
+            self.total_bounds_box = box(*self.total_bounds)
+            
+            if dilated_geometry is not None:
+                self.dilated_bounds = shapely.bounds(dilated_geometry)
+                self.total_dilated_bounds = numpy.array([
+                    numpy.min(self.dilated_bounds[:, 0]),
+                    numpy.min(self.dilated_bounds[:, 1]),
+                    numpy.max(self.dilated_bounds[:, 2]),
+                    numpy.max(self.dilated_bounds[:, 3])
+                ])
+                self.total_dilated_bounds_box = box(*self.total_dilated_bounds)
+            else:
+                self.dilated_bounds = None
+                self.total_dilated_bounds = None
+                self.total_dilated_bounds_box = None
 
     def translate(self, dx: float, dy: float) -> ComponentResult:
         """
         Create a new ComponentResult translated by (dx, dy).
         """
-        # Vectorized translation if possible, else list comp
-        # Shapely 2.x affinity functions still work on single geometries efficiently
+        dxr, dyr = round(dx, 3), round(dy, 3)
+        if (dxr, dyr) == (0.0, 0.0):
+            return self
+        if (dxr, dyr) in self._t_cache:
+            return self._t_cache[(dxr, dyr)]
+
+        # Vectorized translation
         geoms = list(self.geometry)
         num_geom = len(self.geometry)
+        
+        offsets = [num_geom]
         if self.dilated_geometry is not None:
             geoms.extend(self.dilated_geometry)
+        offsets.append(len(geoms))
         
-        from shapely.affinity import translate
-        translated = [translate(p, dx, dy) for p in geoms]
+        if self.proxy_geometry is not None:
+            geoms.extend(self.proxy_geometry)
+        offsets.append(len(geoms))
         
-        new_geom = translated[:num_geom]
-        new_dil = translated[num_geom:] if self.dilated_geometry is not None else None
+        if self.actual_geometry is not None:
+            geoms.extend(self.actual_geometry)
+        offsets.append(len(geoms))
+        
+        import shapely
+        coords = shapely.get_coordinates(geoms)
+        translated = shapely.set_coordinates(geoms, coords + [dx, dy])
+        
+        new_geom = list(translated[:offsets[0]])
+        new_dil = list(translated[offsets[0]:offsets[1]]) if self.dilated_geometry is not None else None
+        new_proxy = list(translated[offsets[1]:offsets[2]]) if self.proxy_geometry is not None else None
+        new_actual = list(translated[offsets[2]:offsets[3]]) if self.actual_geometry is not None else None
         
         new_port = Port(self.end_port.x + dx, self.end_port.y + dy, self.end_port.orientation)
-        return ComponentResult(new_geom, new_port, self.length, new_dil)
-
-
-
-
+        res = ComponentResult(new_geom, new_port, self.length, new_dil, new_proxy, new_actual, skip_bounds=True, move_type=self.move_type)
+        
+        # Optimize: reuse and translate bounds
+        res.bounds = self.bounds + [dx, dy, dx, dy]
+        res.total_bounds = self.total_bounds + [dx, dy, dx, dy]
+        res.total_bounds_box = box(*res.total_bounds)
+        
+        if self.dilated_bounds is not None:
+            res.dilated_bounds = self.dilated_bounds + [dx, dy, dx, dy]
+            res.total_dilated_bounds = self.total_dilated_bounds + [dx, dy, dx, dy]
+            res.total_dilated_bounds_box = box(*res.total_dilated_bounds)
+        else:
+            res.total_dilated_bounds = None
+            res.total_dilated_bounds_box = None
+            
+        self._t_cache[(dxr, dyr)] = res
+        return res
 
 
 class Straight:
@@ -102,19 +144,10 @@ class Straight:
             width: float,
             snap_to_grid: bool = True,
             dilation: float = 0.0,
+            snap_size: float = SEARCH_GRID_SNAP_UM,
             ) -> ComponentResult:
         """
         Generate a straight waveguide segment.
-
-        Args:
-            start_port: Port to start from.
-            length: Requested length.
-            width: Waveguide width.
-            snap_to_grid: Whether to snap the end port to the search grid.
-            dilation: Optional dilation distance for pre-calculating collision geometry.
-
-        Returns:
-            A ComponentResult containing the straight segment.
         """
         rad = numpy.radians(start_port.orientation)
         cos_val = numpy.cos(rad)
@@ -124,8 +157,8 @@ class Straight:
         ey = start_port.y + length * sin_val
 
         if snap_to_grid:
-            ex = snap_search_grid(ex)
-            ey = snap_search_grid(ey)
+            ex = snap_search_grid(ex, snap_size)
+            ey = snap_search_grid(ey, snap_size)
 
         end_port = Port(ex, ey, start_port.orientation)
         actual_length = numpy.sqrt((end_port.x - start_port.x)**2 + (end_port.y - start_port.y)**2)
@@ -159,20 +192,13 @@ class Straight:
             poly_points_dil = (pts_dil @ rot_matrix.T) + [start_port.x, start_port.y]
             dilated_geom = [Polygon(poly_points_dil)]
 
-        return ComponentResult(geometry=geom, end_port=end_port, length=actual_length, dilated_geometry=dilated_geom)
+        # For straight segments, geom IS the actual geometry
+        return ComponentResult(geometry=geom, end_port=end_port, length=actual_length, dilated_geometry=dilated_geom, actual_geometry=geom, move_type='Straight')
 
 
 def _get_num_segments(radius: float, angle_deg: float, sagitta: float = 0.01) -> int:
     """
     Calculate number of segments for an arc to maintain a maximum sagitta.
-
-    Args:
-        radius: Arc radius.
-        angle_deg: Total angle turned.
-        sagitta: Maximum allowed deviation.
-
-    Returns:
-        Minimum number of segments needed.
     """
     if radius <= 0:
         return 1
@@ -196,17 +222,6 @@ def _get_arc_polygons(
         ) -> list[Polygon]:
     """
     Helper to generate arc-shaped polygons using vectorized NumPy operations.
-
-    Args:
-        cx, cy: Center coordinates.
-        radius: Arc radius.
-        width: Waveguide width.
-        t_start, t_end: Start and end angles (radians).
-        sagitta: Geometric fidelity.
-        dilation: Optional dilation to apply directly to the arc.
-
-    Returns:
-        List containing the arc polygon.
     """
     num_segments = _get_num_segments(radius, float(numpy.degrees(abs(t_end - t_start))), sagitta)
     angles = numpy.linspace(t_start, t_end, num_segments + 1)
@@ -227,74 +242,66 @@ def _get_arc_polygons(
 
 
 def _clip_bbox(
-        bbox: Polygon,
         cx: float,
         cy: float,
         radius: float,
         width: float,
-        clip_margin: float,
-        arc_poly: Polygon,
+        t_start: float,
+        t_end: float,
         ) -> Polygon:
     """
-    Clips corners of a bounding box for better collision modeling using direct vertex manipulation.
+    Generates a rotationally invariant bounding polygon for an arc.
     """
-    # Determination of which corners to clip
-    ac = arc_poly.centroid
-    qsx = 1.0 if ac.x >= cx else -1.0
-    qsy = 1.0 if ac.y >= cy else -1.0
-
-    r_out_cut = radius + width / 2.0 + clip_margin
-    r_in_cut = radius - width / 2.0 - clip_margin
-
-    minx, miny, maxx, maxy = bbox.bounds
-    # Initial vertices: [minx,miny], [maxx,miny], [maxx,maxy], [minx,maxy]
-    verts = [
-        numpy.array([minx, miny]),
-        numpy.array([maxx, miny]),
-        numpy.array([maxx, maxy]),
-        numpy.array([minx, maxy])
-    ]
-
-    new_verts = []
-    for p in verts:
-        dx, dy = p[0] - cx, p[1] - cy
-        dist = numpy.sqrt(dx**2 + dy**2)
-
-        # Normal vector components from center to corner
-        sx = 1.0 if dx > 1e-6 else (-1.0 if dx < -1e-6 else qsx)
-        sy = 1.0 if dy > 1e-6 else (-1.0 if dy < -1e-6 else qsy)
-
-        d_line = -1.0
-        if dist > r_out_cut:
-            d_line = r_out_cut * numpy.sqrt(2)
-        elif r_in_cut > 0 and dist < r_in_cut:
-            d_line = r_in_cut
+    sweep = abs(t_end - t_start)
+    if sweep > 2 * numpy.pi:
+        sweep = sweep % (2 * numpy.pi)
+    
+    mid_angle = (t_start + t_end) / 2.0
+    # Handle wrap-around for mid_angle
+    if abs(t_end - t_start) > numpy.pi:
+        mid_angle += numpy.pi
         
-        if d_line > 0:
-            # This corner needs clipping. Replace one vertex with two at intersection of line and edges.
-            # Line: sx*(x-cx) + sy*(y-cy) = d_line
-            # Edge x=px: y = cy + (d_line - sx*(px-cx))/sy
-            # Edge y=py: x = cx + (d_line - sy*(py-cy))/sx
-            try:
-                p_edge_x = numpy.array([p[0], cy + (d_line - sx * (p[0] - cx)) / sy])
-                p_edge_y = numpy.array([cx + (d_line - sy * (p[1] - cy)) / sx, p[1]])
-                # Order matters for polygon winding.
-                # If we are at [minx, miny] and moving CCW towards [maxx, miny]:
-                # If we clip this corner, we should add p_edge_y then p_edge_x (or vice versa depending on orientation)
-                # For simplicity, we can just add both and let Polygon sort it out if it's convex,
-                # but better to be precise.
-                # Since we know the bounding box orientation, we can determine order.
-                # BUT: Difference was safer. Let's try a simpler approach:
-                # Just collect all possible vertices and use convex_hull if it's guaranteed convex.
-                # A clipped bbox is always convex.
-                new_verts.append(p_edge_x)
-                new_verts.append(p_edge_y)
-            except ZeroDivisionError:
-                new_verts.append(p)
-        else:
-            new_verts.append(p)
-
-    return Polygon(new_verts).convex_hull
+    r_out = radius + width / 2.0
+    r_in = max(0.0, radius - width / 2.0)
+    
+    half_sweep = sweep / 2.0
+    
+    # Define vertices in local space (center at 0,0, symmetry axis along +X)
+    # 1. Start Inner
+    # 2. Start Outer
+    # 3. Peak Outer Start (tangent intersection approximation)
+    # 4. Peak Outer End
+    # 5. End Outer
+    # 6. End Inner
+    # 7. Peak Inner (ensures convexity and inner clipping)
+    
+    # To clip the outer corner, we use two peak vertices that follow the arc tighter.
+    cos_hs = numpy.cos(half_sweep)
+    cos_hs2 = numpy.cos(half_sweep / 2.0)
+    tan_hs2 = numpy.tan(half_sweep / 2.0)
+    
+    # Distance to peak from center: r_out / cos(hs/2)
+    # At angles +/- hs/2
+    peak_r = r_out / cos_hs2
+    
+    local_verts = [
+        [r_in * numpy.cos(-half_sweep), r_in * numpy.sin(-half_sweep)],
+        [r_out * numpy.cos(-half_sweep), r_out * numpy.sin(-half_sweep)],
+        [peak_r * numpy.cos(-half_sweep/2), peak_r * numpy.sin(-half_sweep/2)],
+        [peak_r * numpy.cos(half_sweep/2), peak_r * numpy.sin(half_sweep/2)],
+        [r_out * numpy.cos(half_sweep), r_out * numpy.sin(half_sweep)],
+        [r_in * numpy.cos(half_sweep), r_in * numpy.sin(half_sweep)],
+        [r_in, 0.0]
+    ]
+    
+    # Rotate and translate to world space
+    cos_m = numpy.cos(mid_angle)
+    sin_m = numpy.sin(mid_angle)
+    rot = numpy.array([[cos_m, -sin_m], [sin_m, cos_m]])
+    
+    world_verts = (numpy.array(local_verts) @ rot.T) + [cx, cy]
+    
+    return Polygon(world_verts)
 
 
 def _apply_collision_model(
@@ -305,101 +312,123 @@ def _apply_collision_model(
         cx: float = 0.0,
         cy: float = 0.0,
         clip_margin: float = 10.0,
+        t_start: float | None = None,
+        t_end: float | None = None,
         ) -> list[Polygon]:
     """
     Applies the specified collision model to an arc geometry.
-
-    Args:
-        arc_poly: High-fidelity arc.
-        collision_type: Model type or custom polygon.
-        radius: Arc radius.
-        width: Waveguide width.
-        cx, cy: Arc center.
-        clip_margin: Safety margin for clipping.
-
-    Returns:
-        List of polygons representing the collision model.
     """
     if isinstance(collision_type, Polygon):
         return [collision_type]
 
     if collision_type == "arc":
         return [arc_poly]
+    
+    if collision_type == "clipped_bbox" and t_start is not None and t_end is not None:
+        return [_clip_bbox(cx, cy, radius, width, t_start, t_end)]
 
-    # Get bounding box
+    # Bounding box of the high-fidelity arc (fallback for bbox or missing angles)
     minx, miny, maxx, maxy = arc_poly.bounds
-    bbox = box(minx, miny, maxx, maxy)
+    bbox_poly = box(minx, miny, maxx, maxy)
 
     if collision_type == "bbox":
-        return [bbox]
-
-    if collision_type == "clipped_bbox":
-        return [_clip_bbox(bbox, cx, cy, radius, width, clip_margin, arc_poly)]
-
+        return [bbox_poly]
+    
     return [arc_poly]
 
 
 class Bend90:
     """
-    Move generator for 90-degree bends.
+    Move generator for 90-degree waveguide bends.
     """
     @staticmethod
     def generate(
             start_port: Port,
             radius: float,
             width: float,
-            direction: str = "CW",
+            direction: Literal["CW", "CCW"],
             sagitta: float = 0.01,
             collision_type: Literal["arc", "bbox", "clipped_bbox"] | Polygon = "arc",
             clip_margin: float = 10.0,
             dilation: float = 0.0,
+            snap_to_grid: bool = True,
+            snap_size: float = SEARCH_GRID_SNAP_UM,
             ) -> ComponentResult:
         """
         Generate a 90-degree bend.
-
-        Args:
-            start_port: Port to start from.
-            radius: Bend radius.
-            width: Waveguide width.
-            direction: "CW" or "CCW".
-            sagitta: Geometric fidelity.
-            collision_type: Collision model.
-            clip_margin: Margin for clipped_bbox.
-            dilation: Optional dilation distance for pre-calculating collision geometry.
-
-        Returns:
-            A ComponentResult containing the bend.
         """
-        turn_angle = -90 if direction == "CW" else 90
         rad_start = numpy.radians(start_port.orientation)
-        c_angle = rad_start + (numpy.pi / 2 if direction == "CCW" else -numpy.pi / 2)
-        cx = start_port.x + radius * numpy.cos(c_angle)
-        cy = start_port.y + radius * numpy.sin(c_angle)
-        t_start = c_angle + numpy.pi
-        t_end = t_start + (numpy.pi / 2 if direction == "CCW" else -numpy.pi / 2)
-
-        ex = snap_search_grid(cx + radius * numpy.cos(t_end))
-        ey = snap_search_grid(cy + radius * numpy.sin(t_end))
-        end_port = Port(ex, ey, float((start_port.orientation + turn_angle) % 360))
-
-        arc_polys = _get_arc_polygons(cx, cy, radius, width, t_start, t_end, sagitta)
-        collision_polys = _apply_collision_model(
-            arc_polys[0], collision_type, radius, width, cx, cy, clip_margin
-        )
         
+        # Center of the arc
+        if direction == "CCW":
+            cx = start_port.x + radius * numpy.cos(rad_start + numpy.pi / 2)
+            cy = start_port.y + radius * numpy.sin(rad_start + numpy.pi / 2)
+            t_start = rad_start - numpy.pi / 2
+            t_end = t_start + numpy.pi / 2
+            new_ori = (start_port.orientation + 90) % 360
+        else:
+            cx = start_port.x + radius * numpy.cos(rad_start - numpy.pi / 2)
+            cy = start_port.y + radius * numpy.sin(rad_start - numpy.pi / 2)
+            t_start = rad_start + numpy.pi / 2
+            t_end = t_start - numpy.pi / 2
+            new_ori = (start_port.orientation - 90) % 360
+
+        # Snap the end point to the grid
+        ex_raw = cx + radius * numpy.cos(t_end)
+        ey_raw = cy + radius * numpy.sin(t_end)
+        
+        if snap_to_grid:
+            ex = snap_search_grid(ex_raw, snap_size)
+            ey = snap_search_grid(ey_raw, snap_size)
+        else:
+            ex, ey = ex_raw, ey_raw
+        
+        # Slightly adjust radius and t_end to hit snapped point exactly
+        actual_radius = numpy.sqrt((ex - cx)**2 + (ey - cy)**2)
+        
+        t_end_snapped = numpy.arctan2(ey - cy, ex - cx)
+        # Ensure directionality and approx 90 degree sweep
+        if direction == "CCW":
+            while t_end_snapped <= t_start:
+                t_end_snapped += 2 * numpy.pi
+            while t_end_snapped > t_start + numpy.pi:
+                t_end_snapped -= 2 * numpy.pi
+        else:
+            while t_end_snapped >= t_start:
+                t_end_snapped -= 2 * numpy.pi
+            while t_end_snapped < t_start - numpy.pi:
+                t_end_snapped += 2 * numpy.pi
+        t_end = t_end_snapped
+
+        end_port = Port(ex, ey, new_ori)
+
+        arc_polys = _get_arc_polygons(cx, cy, actual_radius, width, t_start, t_end, sagitta)
+        collision_polys = _apply_collision_model(
+            arc_polys[0], collision_type, actual_radius, width, cx, cy, clip_margin, t_start, t_end
+        )
+
+        proxy_geom = None
+        if collision_type == "arc":
+            # Auto-generate a clipped_bbox proxy for tiered collision checks
+            proxy_geom = _apply_collision_model(
+                arc_polys[0], "clipped_bbox", actual_radius, width, cx, cy, clip_margin, t_start, t_end
+            )
+
         dilated_geom = None
         if dilation > 0:
             if collision_type == "arc":
-                dilated_geom = _get_arc_polygons(cx, cy, radius, width, t_start, t_end, sagitta, dilation=dilation)
+                dilated_geom = _get_arc_polygons(cx, cy, actual_radius, width, t_start, t_end, sagitta, dilation=dilation)
             else:
-                # For bbox or clipped_bbox, buffer the model itself (which is simpler than buffering the high-fidelity arc)
                 dilated_geom = [p.buffer(dilation) for p in collision_polys]
 
         return ComponentResult(
             geometry=collision_polys, 
             end_port=end_port, 
-            length=radius * numpy.pi / 2.0,
-            dilated_geometry=dilated_geom
+            length=actual_radius * numpy.abs(t_end - t_start),
+            dilated_geometry=dilated_geom,
+            proxy_geometry=proxy_geom,
+            actual_geometry=arc_polys,
+            move_type='Bend90'
         )
 
 
@@ -417,68 +446,90 @@ class SBend:
             collision_type: Literal["arc", "bbox", "clipped_bbox"] | Polygon = "arc",
             clip_margin: float = 10.0,
             dilation: float = 0.0,
+            snap_to_grid: bool = True,
+            snap_size: float = SEARCH_GRID_SNAP_UM,
             ) -> ComponentResult:
         """
         Generate a parametric S-bend (two tangent arcs).
-
-        Args:
-            start_port: Port to start from.
-            offset: Lateral offset.
-            radius: Arc radii.
-            width: Waveguide width.
-            sagitta: Geometric fidelity.
-            collision_type: Collision model.
-            clip_margin: Margin for clipped_bbox.
-            dilation: Optional dilation distance for pre-calculating collision geometry.
-
-        Returns:
-            A ComponentResult containing the S-bend.
         """
         if abs(offset) >= 2 * radius:
             raise ValueError(f"SBend offset {offset} must be less than 2*radius {2 * radius}")
 
-        theta = numpy.arccos(1 - abs(offset) / (2 * radius))
-        dx = 2 * radius * numpy.sin(theta)
-        dy = offset
+        theta_init = numpy.arccos(1 - abs(offset) / (2 * radius))
+        dx_init = 2 * radius * numpy.sin(theta_init)
         rad_start = numpy.radians(start_port.orientation)
-        ex = snap_search_grid(start_port.x + dx * numpy.cos(rad_start) - dy * numpy.sin(rad_start))
-        ey = snap_search_grid(start_port.y + dx * numpy.sin(rad_start) + dy * numpy.cos(rad_start))
+        
+        # Target point
+        ex_raw = start_port.x + dx_init * numpy.cos(rad_start) - offset * numpy.sin(rad_start)
+        ey_raw = start_port.y + dx_init * numpy.sin(rad_start) + offset * numpy.cos(rad_start)
+
+        if snap_to_grid:
+            ex = snap_search_grid(ex_raw, snap_size)
+            ey = snap_search_grid(ey_raw, snap_size)
+        else:
+            ex, ey = ex_raw, ey_raw
+
         end_port = Port(ex, ey, start_port.orientation)
 
-        direction = 1 if offset > 0 else -1
+        # Solve for theta and radius that hit (ex, ey) exactly
+        local_dx = (ex - start_port.x) * numpy.cos(rad_start) + (ey - start_port.y) * numpy.sin(rad_start)
+        local_dy = -(ex - start_port.x) * numpy.sin(rad_start) + (ey - start_port.y) * numpy.cos(rad_start)
+        
+        # tan(theta / 2) = local_dy / local_dx
+        theta = 2 * numpy.arctan2(abs(local_dy), local_dx)
+        
+        if abs(theta) < 1e-9:
+             # Practically straight, but offset implies we need a bend.
+             # If offset is also tiny, return a straight?
+             if abs(offset) < 1e-6:
+                 # Degenerate case: effectively straight
+                 return Straight.generate(start_port, numpy.sqrt(local_dx**2 + local_dy**2), width, snap_to_grid=False, dilation=dilation)
+             raise ValueError("SBend calculation failed: theta close to zero")
+
+        # Avoid division by zero if theta is 0 (though unlikely due to offset check)
+        denom = (2 * (1 - numpy.cos(theta)))
+        if abs(denom) < 1e-9:
+             raise ValueError("SBend calculation failed: radius denominator zero")
+             
+        actual_radius = abs(local_dy) / denom
+        
+        # Limit radius to prevent giant arcs
+        if actual_radius > 100000.0:
+             raise ValueError("SBend calculation failed: radius too large")
+
+        direction = 1 if local_dy > 0 else -1
         c1_angle = rad_start + direction * numpy.pi / 2
-        cx1 = start_port.x + radius * numpy.cos(c1_angle)
-        cy1 = start_port.y + radius * numpy.sin(c1_angle)
+        cx1 = start_port.x + actual_radius * numpy.cos(c1_angle)
+        cy1 = start_port.y + actual_radius * numpy.sin(c1_angle)
         ts1, te1 = c1_angle + numpy.pi, c1_angle + numpy.pi + direction * theta
 
-        ex_raw = start_port.x + dx * numpy.cos(rad_start) - dy * numpy.sin(rad_start)
-        ey_raw = start_port.y + dx * numpy.sin(rad_start) + dy * numpy.cos(rad_start)
         c2_angle = rad_start - direction * numpy.pi / 2
-        cx2 = ex_raw + radius * numpy.cos(c2_angle)
-        cy2 = ey_raw + radius * numpy.sin(c2_angle)
+        cx2 = ex + actual_radius * numpy.cos(c2_angle)
+        cy2 = ey + actual_radius * numpy.sin(c2_angle)
         te2 = c2_angle + numpy.pi
         ts2 = te2 + direction * theta
 
-        arc1 = _get_arc_polygons(cx1, cy1, radius, width, ts1, te1, sagitta)[0]
-        arc2 = _get_arc_polygons(cx2, cy2, radius, width, ts2, te2, sagitta)[0]
+        arc1 = _get_arc_polygons(cx1, cy1, actual_radius, width, ts1, te1, sagitta)[0]
+        arc2 = _get_arc_polygons(cx2, cy2, actual_radius, width, ts2, te2, sagitta)[0]
+        arc_polys = [arc1, arc2]
 
-        if collision_type == "clipped_bbox":
-            col1 = _apply_collision_model(arc1, collision_type, radius, width, cx1, cy1, clip_margin)[0]
-            col2 = _apply_collision_model(arc2, collision_type, radius, width, cx2, cy2, clip_margin)[0]
-            # Optimization: keep as list instead of unary_union for search efficiency
-            collision_polys = [col1, col2]
-        else:
-            # For other models, we can either combine or keep separate. 
-            # Keeping separate is generally better for CollisionEngine.
-            col1 = _apply_collision_model(arc1, collision_type, radius, width, cx1, cy1, clip_margin)[0]
-            col2 = _apply_collision_model(arc2, collision_type, radius, width, cx2, cy2, clip_margin)[0]
-            collision_polys = [col1, col2]
+        # Use the provided collision model for primary geometry
+        col1 = _apply_collision_model(arc1, collision_type, actual_radius, width, cx1, cy1, clip_margin, ts1, te1)[0]
+        col2 = _apply_collision_model(arc2, collision_type, actual_radius, width, cx2, cy2, clip_margin, ts2, te2)[0]
+        collision_polys = [col1, col2]
+
+        proxy_geom = None
+        if collision_type == "arc":
+            # Auto-generate proxies
+            p1 = _apply_collision_model(arc1, "clipped_bbox", actual_radius, width, cx1, cy1, clip_margin, ts1, te1)[0]
+            p2 = _apply_collision_model(arc2, "clipped_bbox", actual_radius, width, cx2, cy2, clip_margin, ts2, te2)[0]
+            proxy_geom = [p1, p2]
 
         dilated_geom = None
         if dilation > 0:
             if collision_type == "arc":
-                d1 = _get_arc_polygons(cx1, cy1, radius, width, ts1, te1, sagitta, dilation=dilation)[0]
-                d2 = _get_arc_polygons(cx2, cy2, radius, width, ts2, te2, sagitta, dilation=dilation)[0]
+                d1 = _get_arc_polygons(cx1, cy1, actual_radius, width, ts1, te1, sagitta, dilation=dilation)[0]
+                d2 = _get_arc_polygons(cx2, cy2, actual_radius, width, ts2, te2, sagitta, dilation=dilation)[0]
                 dilated_geom = [d1, d2]
             else:
                 dilated_geom = [p.buffer(dilation) for p in collision_polys]
@@ -486,6 +537,9 @@ class SBend:
         return ComponentResult(
             geometry=collision_polys, 
             end_port=end_port, 
-            length=2 * radius * theta,
-            dilated_geometry=dilated_geom
+            length=2 * actual_radius * theta,
+            dilated_geometry=dilated_geom,
+            proxy_geometry=proxy_geom,
+            actual_geometry=arc_polys,
+            move_type='SBend'
         )

inire/geometry/primitives.py

@@ -10,14 +10,8 @@ GRID_SNAP_UM = 0.001
 def snap_nm(value: float) -> float:
     """
     Snap a coordinate to the nearest 1nm (0.001 um).
-
-    Args:
-        value: Coordinate value to snap.
-
-    Returns:
-        Snapped coordinate value.
     """
-    return round(value / GRID_SNAP_UM) * GRID_SNAP_UM
+    return round(value * 1000) / 1000
 
 
 class Port:
@@ -26,39 +20,15 @@ class Port:
     """
     __slots__ = ('x', 'y', 'orientation')
 
-    x: float
-    """ x-coordinate in micrometers """
-
-    y: float
-    """ y-coordinate in micrometers """
-
-    orientation: float
-    """ Orientation in degrees: 0, 90, 180, 270 """
-
     def __init__(
             self,
             x: float,
             y: float,
             orientation: float,
             ) -> None:
-        """
-        Initialize and snap a Port.
-
-        Args:
-            x: Initial x-coordinate.
-            y: Initial y-coordinate.
-            orientation: Initial orientation in degrees.
-        """
-        # Snap x, y to 1nm
         self.x = snap_nm(x)
         self.y = snap_nm(y)
-
-        # Ensure orientation is one of {0, 90, 180, 270}
-        norm_orientation = int(round(orientation)) % 360
-        if norm_orientation not in {0, 90, 180, 270}:
-            norm_orientation = (round(norm_orientation / 90) * 90) % 360
-
-        self.orientation = float(norm_orientation)
+        self.orientation = float(orientation % 360)
 
     def __repr__(self) -> str:
         return f'Port(x={self.x}, y={self.y}, orientation={self.orientation})'
@@ -77,14 +47,6 @@ class Port:
 def translate_port(port: Port, dx: float, dy: float) -> Port:
     """
     Translate a port by (dx, dy).
-
-    Args:
-        port: Port to translate.
-        dx: x-offset.
-        dy: y-offset.
-
-    Returns:
-        A new translated Port.
     """
     return Port(port.x + dx, port.y + dy, port.orientation)
 
@@ -92,20 +54,12 @@ def translate_port(port: Port, dx: float, dy: float) -> Port:
 def rotate_port(port: Port, angle: float, origin: tuple[float, float] = (0, 0)) -> Port:
     """
     Rotate a port by a multiple of 90 degrees around an origin.
-
-    Args:
-        port: Port to rotate.
-        angle: Angle to rotate by (degrees).
-        origin: (x, y) origin to rotate around.
-
-    Returns:
-        A new rotated Port.
     """
     ox, oy = origin
     px, py = port.x, port.y
 
     rad = numpy.radians(angle)
-    qx = ox + numpy.cos(rad) * (px - ox) - numpy.sin(rad) * (py - oy)
-    qy = oy + numpy.sin(rad) * (px - ox) + numpy.cos(rad) * (py - oy)
+    qx = snap_nm(ox + numpy.cos(rad) * (px - ox) - numpy.sin(rad) * (py - oy))
+    qy = snap_nm(oy + numpy.sin(rad) * (px - ox) + numpy.cos(rad) * (py - oy))
 
     return Port(qx, qy, port.orientation + angle)

inire/router/astar.py

@@ -7,10 +7,12 @@ from typing import TYPE_CHECKING, Literal, Any
 import rtree
 
 import numpy
+import shapely
 
-from inire.geometry.components import Bend90, SBend, Straight
+from inire.geometry.components import Bend90, SBend, Straight, SEARCH_GRID_SNAP_UM, snap_search_grid
 from inire.geometry.primitives import Port
 from inire.router.config import RouterConfig
+from inire.router.visibility import VisibilityManager
 
 if TYPE_CHECKING:
     from inire.geometry.components import ComponentResult
@@ -19,38 +21,11 @@ if TYPE_CHECKING:
 logger = logging.getLogger(__name__)
 
 
-@functools.total_ordering
 class AStarNode:
     """
-    A node in the A* search graph.
+    A node in the A* search tree.
     """
-    __slots__ = ('port', 'g_cost', 'h_cost', 'f_cost', 'parent', 'component_result', 'count', 'path_bbox')
-
-    port: Port
-    """ Port representing the state at this node """
-
-    g_cost: float
-    """ Actual cost from start to this node """
-
-    h_cost: float
-    """ Heuristic cost from this node to target """
-
-    f_cost: float
-    """ Total estimated cost (g + h) """
-
-    parent: AStarNode | None
-    """ Parent node in the search tree """
-
-    component_result: ComponentResult | None
-    """ The component move that led to this node """
-
-    count: int
-    """ Unique insertion order for tie-breaking """
-
-    path_bbox: tuple[float, float, float, float] | None
-    """ Bounding box of the entire path up to this node """
-
-    _count = 0
+    __slots__ = ('port', 'g_cost', 'h_cost', 'f_cost', 'parent', 'component_result')
 
     def __init__(
             self,
@@ -66,127 +41,71 @@ class AStarNode:
         self.f_cost = g_cost + h_cost
         self.parent = parent
         self.component_result = component_result
-        self.count = AStarNode._count
-        AStarNode._count += 1
-
-        # Calculate path_bbox
-        if parent is None:
-            self.path_bbox = None
-        else:
-            # Union of parent's bbox and current move's bbox
-            if component_result:
-                # Merge all polygon bounds in the result
-                minx, miny, maxx, maxy = 1e15, 1e15, -1e15, -1e15
-                for b in component_result.dilated_bounds if component_result.dilated_bounds is not None else component_result.bounds:
-                    minx = min(minx, b[0])
-                    miny = min(miny, b[1])
-                    maxx = max(maxx, b[2])
-                    maxy = max(maxy, b[3])
-                
-                if parent.path_bbox:
-                    self.path_bbox = (
-                        min(minx, parent.path_bbox[0]),
-                        min(miny, parent.path_bbox[1]),
-                        max(maxx, parent.path_bbox[2]),
-                        max(maxy, parent.path_bbox[3])
-                    )
-                else:
-                    self.path_bbox = (minx, miny, maxx, maxy)
-            else:
-                self.path_bbox = parent.path_bbox
-
-
-
 
     def __lt__(self, other: AStarNode) -> bool:
-        # Tie-breaking: lower f first, then lower h, then order
-        if abs(self.f_cost - other.f_cost) > 1e-9:
-            return self.f_cost < other.f_cost
-        if abs(self.h_cost - other.h_cost) > 1e-9:
-            return self.h_cost < other.h_cost
-        return self.count < other.count
-
-    def __eq__(self, other: object) -> bool:
-        if not isinstance(other, AStarNode):
+        if self.f_cost < other.f_cost - 1e-6:
+            return True
+        if self.f_cost > other.f_cost + 1e-6:
             return False
-        return self.count == other.count
+        return self.h_cost < other.h_cost
 
 
 class AStarRouter:
     """
-    Hybrid State-Lattice A* Router.
+    Waveguide router based on sparse A* search.
     """
-    __slots__ = ('cost_evaluator', 'config', 'node_limit', 'total_nodes_expanded', '_collision_cache', '_move_cache', '_self_dilation')
-
-    cost_evaluator: CostEvaluator
-    """ The evaluator for path and proximity costs """
-
-    config: RouterConfig
-    """ Search configuration parameters """
-
-    node_limit: int
-    """ Maximum nodes to expand before failure """
-
-    total_nodes_expanded: int
-    """ Counter for debugging/profiling """
-
-    _collision_cache: dict[tuple[float, float, float, str, float, str], bool]
-    """ Internal cache for move collision checks """
-
-    _move_cache: dict[tuple[Any, ...], ComponentResult]
-    """ Internal cache for component generation """
-
-    _self_dilation: float
-    """ Cached dilation value for collision checks (clearance / 2.0) """
-
-    def __init__(
-            self,
-            cost_evaluator: CostEvaluator,
-            node_limit: int = 1000000,
-            straight_lengths: list[float] | None = None,
-            bend_radii: list[float] | None = None,
-            sbend_offsets: list[float] | None = None,
-            sbend_radii: list[float] | None = None,
-            snap_to_target_dist: float = 20.0,
-            bend_penalty: float = 50.0,
-            sbend_penalty: float = 100.0,
-            bend_collision_type: Literal['arc', 'bbox', 'clipped_bbox'] = 'arc',
-            bend_clip_margin: float = 10.0,
-            ) -> None:
-        """
-        Initialize the A* Router.
-
-        Args:
-            cost_evaluator: Path cost evaluator.
-            node_limit: Node expansion limit.
-            straight_lengths: Allowed straight lengths (um).
-            bend_radii: Allowed 90-deg radii (um).
-            sbend_offsets: Allowed S-bend lateral offsets (um).
-            sbend_radii: Allowed S-bend radii (um).
-            snap_to_target_dist: Radius for target lookahead (um).
-            bend_penalty: Penalty for 90-degree turns.
-            sbend_penalty: Penalty for S-bends.
-            bend_collision_type: Collision model for bends.
-            bend_clip_margin: Margin for clipped_bbox model.
-        """
+    def __init__(self, cost_evaluator: CostEvaluator, node_limit: int | None = None, **kwargs) -> None:
         self.cost_evaluator = cost_evaluator
-        self.config = RouterConfig(
-            node_limit=node_limit,
-            straight_lengths=straight_lengths if straight_lengths is not None else [1.0, 5.0, 25.0, 100.0],
-            bend_radii=bend_radii if bend_radii is not None else [10.0],
-            sbend_offsets=sbend_offsets if sbend_offsets is not None else [-5.0, -2.0, 2.0, 5.0],
-            sbend_radii=sbend_radii if sbend_radii is not None else [10.0],
-            snap_to_target_dist=snap_to_target_dist,
-            bend_penalty=bend_penalty,
-            sbend_penalty=sbend_penalty,
-            bend_collision_type=bend_collision_type,
-            bend_clip_margin=bend_clip_margin,
-        )
+        self.config = RouterConfig()
+        
+        if node_limit is not None:
+            self.config.node_limit = node_limit
+            
+        for k, v in kwargs.items():
+            if hasattr(self.config, k):
+                setattr(self.config, k, v)
+        
         self.node_limit = self.config.node_limit
+        
+        # Visibility Manager for sparse jumps
+        self.visibility_manager = VisibilityManager(self.cost_evaluator.collision_engine)
+        
+        self._hard_collision_set: set[tuple] = set()
+        self._congestion_cache: dict[tuple, int] = {}
+        self._static_safe_cache: set[tuple] = set()
+        self._move_cache: dict[tuple, ComponentResult] = {}
+
         self.total_nodes_expanded = 0
-        self._collision_cache = {}
-        self._move_cache = {}
-        self._self_dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+        self.last_expanded_nodes: list[tuple[float, float, float]] = []
+        
+        self.metrics = {
+            'nodes_expanded': 0,
+            'moves_generated': 0,
+            'moves_added': 0,
+            'pruned_closed_set': 0,
+            'pruned_hard_collision': 0,
+            'pruned_cost': 0
+        }
+
+    def reset_metrics(self) -> None:
+        """ Reset all performance counters. """
+        for k in self.metrics:
+            self.metrics[k] = 0
+        self.cost_evaluator.collision_engine.reset_metrics()
+
+    def get_metrics_summary(self) -> str:
+        """ Return a human-readable summary of search performance. """
+        m = self.metrics
+        c = self.cost_evaluator.collision_engine.get_metrics_summary()
+        return (f"Search Performance: \n"
+                f"  Nodes Expanded: {m['nodes_expanded']}\n"
+                f"  Moves: Generated={m['moves_generated']}, Added={m['moves_added']}\n"
+                f"  Pruning: ClosedSet={m['pruned_closed_set']}, HardColl={m['pruned_hard_collision']}, Cost={m['pruned_cost']}\n"
+                f"  {c}")
+
+    @property
+    def _self_dilation(self) -> float:
+        return self.cost_evaluator.collision_engine.clearance / 2.0
 
     def route(
             self,
@@ -195,51 +114,54 @@ class AStarRouter:
             net_width: float,
             net_id: str = 'default',
             bend_collision_type: Literal['arc', 'bbox', 'clipped_bbox'] | None = None,
+            return_partial: bool = False,
+            store_expanded: bool = False,
+            skip_congestion: bool = False,
             ) -> list[ComponentResult] | None:
         """
         Route a single net using A*.
-
-        Args:
-            start: Starting port.
-            target: Target port.
-            net_width: Waveguide width (um).
-            net_id: Optional net identifier.
-            bend_collision_type: Override collision model for this route.
-
-        Returns:
-            List of moves forming the path, or None if failed.
         """
+        self._congestion_cache.clear()
+        if store_expanded:
+            self.last_expanded_nodes = []
+
         if bend_collision_type is not None:
             self.config.bend_collision_type = bend_collision_type
             
-        self._collision_cache.clear()
         open_set: list[AStarNode] = []
-        # Key: (x, y, orientation) rounded to 1nm
-        closed_set: set[tuple[float, float, float]] = set()
+        snap = self.config.snap_size
+        
+        # (x_grid, y_grid, orientation_grid) -> min_g_cost
+        closed_set: dict[tuple[int, int, int], float] = {}
 
         start_node = AStarNode(start, 0.0, self.cost_evaluator.h_manhattan(start, target))
         heapq.heappush(open_set, start_node)
 
+        best_node = start_node
         nodes_expanded = 0
+        
+        node_limit = self.node_limit
 
         while open_set:
-            if nodes_expanded >= self.node_limit:
-                logger.warning(f'  AStar failed: node limit {self.node_limit} reached.')
-                return None
+            if nodes_expanded >= node_limit:
+                return self._reconstruct_path(best_node) if return_partial else None
 
             current = heapq.heappop(open_set)
+            
+            if current.h_cost < best_node.h_cost:
+                best_node = current
 
-            # Prune if already visited
-            state = (round(current.port.x, 3), round(current.port.y, 3), round(current.port.orientation, 2))
-            if state in closed_set:
+            state = (int(round(current.port.x / snap)), int(round(current.port.y / snap)), int(round(current.port.orientation / 1.0)))
+            if state in closed_set and closed_set[state] <= current.g_cost + 1e-6:
                 continue
-            closed_set.add(state)
+            closed_set[state] = current.g_cost
+
+            if store_expanded:
+                self.last_expanded_nodes.append((current.port.x, current.port.y, current.port.orientation))
 
             nodes_expanded += 1
             self.total_nodes_expanded += 1
-
-            if nodes_expanded % 5000 == 0:
-                logger.info(f'Nodes expanded: {nodes_expanded}, current: {current.port}, g: {current.g_cost:.1f}')
+            self.metrics['nodes_expanded'] += 1
 
             # Check if we reached the target exactly
             if (abs(current.port.x - target.x) < 1e-6 and
@@ -248,9 +170,9 @@ class AStarRouter:
                 return self._reconstruct_path(current)
 
             # Expansion
-            self._expand_moves(current, target, net_width, net_id, open_set, closed_set)
+            self._expand_moves(current, target, net_width, net_id, open_set, closed_set, snap, nodes_expanded, skip_congestion=skip_congestion)
 
-        return None
+        return self._reconstruct_path(best_node) if return_partial else None
 
     def _expand_moves(
             self,
@@ -259,148 +181,210 @@ class AStarRouter:
             net_width: float,
             net_id: str,
             open_set: list[AStarNode],
-            closed_set: set[tuple[float, float, float]],
+            closed_set: dict[tuple[int, int, int], float],
+            snap: float = 1.0,
+            nodes_expanded: int = 0,
+            skip_congestion: bool = False,
             ) -> None:
-        # 1. Snap-to-Target Look-ahead
-        dist = numpy.sqrt((current.port.x - target.x)**2 + (current.port.y - target.y)**2)
-        if dist < self.config.snap_to_target_dist:
-            # A. Try straight exact reach
-            if abs(current.port.orientation - target.orientation) < 0.1:
-                rad = numpy.radians(current.port.orientation)
-                dx = target.x - current.port.x
-                dy = target.y - current.port.y
-                proj = dx * numpy.cos(rad) + dy * numpy.sin(rad)
-                perp = -dx * numpy.sin(rad) + dy * numpy.cos(rad)
-                if proj > 0 and abs(perp) < 1e-6:
-                    res = Straight.generate(current.port, proj, net_width, snap_to_grid=False, dilation=self._self_dilation)
-                    self._add_node(current, res, target, net_width, net_id, open_set, closed_set, 'SnapStraight')
-
-            # B. Try SBend exact reach
-            if abs(current.port.orientation - target.orientation) < 0.1:
-                rad = numpy.radians(current.port.orientation)
-                dx = target.x - current.port.x
-                dy = target.y - current.port.y
-                proj = dx * numpy.cos(rad) + dy * numpy.sin(rad)
-                perp = -dx * numpy.sin(rad) + dy * numpy.cos(rad)
-                if proj > 0 and 0.5 <= abs(perp) < 20.0:
-                    for radius in self.config.sbend_radii:
-                        try:
-                            res = SBend.generate(
-                                current.port,
-                                perp,
-                                radius,
-                                net_width,
-                                collision_type=self.config.bend_collision_type,
-                                clip_margin=self.config.bend_clip_margin,
-                                dilation=self._self_dilation
-                            )
-                            self._add_node(current, res, target, net_width, net_id, open_set, closed_set, 'SnapSBend', move_radius=radius)
-                        except ValueError:
-                            pass
-
-        # Move Cache
         cp = current.port
-        base_ori = round(cp.orientation % 360, 2)
-        state_key = (round(cp.x, 3), round(cp.y, 3), base_ori)
+        base_ori = round(cp.orientation, 2)
+        dx_t = target.x - cp.x
+        dy_t = target.y - cp.y
+        dist_sq = dx_t*dx_t + dy_t*dy_t
+        
+        rad = numpy.radians(base_ori)
+        cos_v, sin_v = numpy.cos(rad), numpy.sin(rad)
+        # 1. DIRECT JUMP TO TARGET (Priority 1)
+        proj_t = dx_t * cos_v + dy_t * sin_v
+        perp_t = -dx_t * sin_v + dy_t * cos_v
 
-        # 2. Lattice Straights
-        lengths = self.config.straight_lengths
-        if dist < 5.0:
-             fine_steps = [0.1, 0.5]
-             lengths = sorted(set(lengths + fine_steps))
+        # A. Straight Jump
+        if proj_t > 0 and abs(perp_t) < 1e-3 and abs(cp.orientation - target.orientation) < 0.1:
+             max_reach = self.cost_evaluator.collision_engine.ray_cast(cp, base_ori, proj_t + 1.0)
+             if max_reach >= proj_t - 0.01:
+                  self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{proj_t}', 'S', (proj_t,), skip_congestion, skip_static=True, snap_to_grid=False)
 
-        for length in lengths:
-            # Level 1: Absolute cache (exact location)
-            abs_key = (state_key, 'S', length, net_width)
-            if abs_key in self._move_cache:
-                res = self._move_cache[abs_key]
-            else:
-                # Level 2: Relative cache (orientation only)
-                # Dilation is now 0.0 for caching to save translation time.
-                # It will be calculated lazily in _add_node if needed.
-                rel_key = (base_ori, 'S', length, net_width, 0.0)
-                if rel_key in self._move_cache:
-                    res_rel = self._move_cache[rel_key]
-                    # Check closed set before translating
-                    ex = res_rel.end_port.x + cp.x
-                    ey = res_rel.end_port.y + cp.y
-                    end_state = (round(ex, 3), round(ey, 3), round(res_rel.end_port.orientation, 2))
-                    if end_state in closed_set:
-                        continue
-                    res = res_rel.translate(cp.x, cp.y)
-                else:
-                    res_rel = Straight.generate(Port(0, 0, base_ori), length, net_width, dilation=0.0)
-                    self._move_cache[rel_key] = res_rel
-                    res = res_rel.translate(cp.x, cp.y)
-                self._move_cache[abs_key] = res
-            self._add_node(current, res, target, net_width, net_id, open_set, closed_set, f'S{length}')
+        # B. SBend Jump (Direct to Target)
+        if self.config.use_analytical_sbends and proj_t > 0 and abs(cp.orientation - target.orientation) < 0.1 and abs(perp_t) > 1e-3:
+             # Calculate required radius to hit target exactly: R = (dx^2 + dy^2) / (4*|dy|)
+             req_radius = (proj_t**2 + perp_t**2) / (4.0 * abs(perp_t))
+             
+             min_radius = min(self.config.sbend_radii) if self.config.sbend_radii else 50.0
+             
+             if req_radius >= min_radius:
+                  # We can hit it exactly!
+                  self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'SB_Direct_R{req_radius:.1f}', 'SB', (perp_t, req_radius), skip_congestion, snap_to_grid=False)
+             else:
+                  # Required radius is too small. We must use a larger radius and some straight segments.
+                  # A* will handle this through Priority 3 SBends + Priority 2 Straights.
+                  pass
+
+                  # In super sparse mode, we can return here, but A* needs other options for optimality.
+                  # return 
+
+        # 2. VISIBILITY JUMPS & MAX REACH (Priority 2)
+        max_reach = self.cost_evaluator.collision_engine.ray_cast(cp, base_ori, self.config.max_straight_length)
+        
+        straight_lengths = set()
+        if max_reach > self.config.min_straight_length:
+             # milestone 1: exactly at max_reach (touching)
+             straight_lengths.add(snap_search_grid(max_reach, snap))
+             # milestone 2: space to turn before collision
+             for radius in self.config.bend_radii:
+                  if max_reach > radius + self.config.min_straight_length:
+                       straight_lengths.add(snap_search_grid(max_reach - radius, snap))
+             
+             # milestone 3: small buffer for tight maneuvering
+             if max_reach > self.config.min_straight_length + 5.0:
+                  straight_lengths.add(snap_search_grid(max_reach - 5.0, snap))
+
+        visible_corners = self.visibility_manager.get_visible_corners(cp, max_dist=max_reach)
+        for cx, cy, dist in visible_corners:
+            proj = (cx - cp.x) * cos_v + (cy - cp.y) * sin_v
+            if proj > self.config.min_straight_length:
+                 straight_lengths.add(snap_search_grid(proj, snap))
+        
+        # ALWAYS include the min length for maneuvering
+        straight_lengths.add(self.config.min_straight_length)
+        
+        # If the jump is long, add an intermediate point to allow more flexible turning
+        if max_reach > self.config.min_straight_length * 4:
+             straight_lengths.add(snap_search_grid(max_reach / 2.0, snap))
+
+        # Target alignment logic (for turning towards target) - Keep this as it's high value
+        if abs(base_ori % 180) < 0.1: # Horizontal
+             target_dist = abs(target.x - cp.x)
+             if target_dist <= max_reach and target_dist > self.config.min_straight_length:
+                  straight_lengths.add(snap_search_grid(target_dist, snap))
+             
+             # Space for turning: target_dist - R and target_dist - 2R
+             for radius in self.config.bend_radii:
+                  l1 = target_dist - radius
+                  if l1 > self.config.min_straight_length:
+                       s_l1 = snap_search_grid(l1, snap)
+                       if s_l1 <= max_reach and s_l1 > 0.1:
+                            straight_lengths.add(s_l1)
+                  
+                  l2 = target_dist - 2 * radius
+                  if l2 > self.config.min_straight_length:
+                       s_l2 = snap_search_grid(l2, snap)
+                       if s_l2 <= max_reach and s_l2 > 0.1:
+                            straight_lengths.add(s_l2)
+        else: # Vertical
+             target_dist = abs(target.y - cp.y)
+             if target_dist <= max_reach and target_dist > self.config.min_straight_length:
+                  straight_lengths.add(snap_search_grid(target_dist, snap))
+             
+             # Space for turning: target_dist - R and target_dist - 2R
+             for radius in self.config.bend_radii:
+                  l1 = target_dist - radius
+                  if l1 > self.config.min_straight_length:
+                       s_l1 = snap_search_grid(l1, snap)
+                       if s_l1 <= max_reach and s_l1 > 0.1:
+                            straight_lengths.add(s_l1)
+                  
+                  l2 = target_dist - 2 * radius
+                  if l2 > self.config.min_straight_length:
+                       s_l2 = snap_search_grid(l2, snap)
+                       if s_l2 <= max_reach and s_l2 > 0.1:
+                            straight_lengths.add(s_l2)
+
+        # NO standard samples here! Only milestones.
+
+        for length in sorted(straight_lengths, reverse=True):
+            # Trust ray_cast: these lengths are <= max_reach
+            self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'S{length}', 'S', (length,), skip_congestion, skip_static=True)
+
+        # 3. BENDS & SBENDS (Priority 3)
+        angle_to_target = numpy.degrees(numpy.arctan2(target.y - cp.y, target.x - cp.x))
+        allow_backwards = (dist_sq < 150*150)
 
-        # 3. Lattice Bends
         for radius in self.config.bend_radii:
             for direction in ['CW', 'CCW']:
-                abs_key = (state_key, 'B', radius, direction, net_width, self.config.bend_collision_type)
-                if abs_key in self._move_cache:
-                    res = self._move_cache[abs_key]
-                else:
-                    rel_key = (base_ori, 'B', radius, direction, net_width, self.config.bend_collision_type, 0.0)
-                    if rel_key in self._move_cache:
-                        res_rel = self._move_cache[rel_key]
-                        # Check closed set before translating
-                        ex = res_rel.end_port.x + cp.x
-                        ey = res_rel.end_port.y + cp.y
-                        end_state = (round(ex, 3), round(ey, 3), round(res_rel.end_port.orientation, 2))
-                        if end_state in closed_set:
-                            continue
-                        res = res_rel.translate(cp.x, cp.y)
-                    else:
-                        res_rel = Bend90.generate(
-                            Port(0, 0, base_ori),
-                            radius,
-                            net_width,
-                            direction,
-                            collision_type=self.config.bend_collision_type,
-                            clip_margin=self.config.bend_clip_margin,
-                            dilation=0.0
-                        )
-                        self._move_cache[rel_key] = res_rel
-                        res = res_rel.translate(cp.x, cp.y)
-                    self._move_cache[abs_key] = res
-                self._add_node(current, res, target, net_width, net_id, open_set, closed_set, f'B{radius}{direction}', move_radius=radius)
+                if not allow_backwards:
+                    turn = 90 if direction == 'CCW' else -90
+                    new_ori = (cp.orientation + turn) % 360
+                    new_diff = (angle_to_target - new_ori + 180) % 360 - 180
+                    if abs(new_diff) > 135:
+                        continue
+                self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'B{radius}{direction}', 'B', (radius, direction), skip_congestion)
 
-        # 4. Discrete SBends
-        for offset in self.config.sbend_offsets:
-            for radius in self.config.sbend_radii:
-                abs_key = (state_key, 'SB', offset, radius, net_width, self.config.bend_collision_type)
-                if abs_key in self._move_cache:
-                    res = self._move_cache[abs_key]
+        if dist_sq < 400*400:
+            offsets = set(self.config.sbend_offsets)
+            dx_local = (target.x - cp.x) * cos_v + (target.y - cp.y) * sin_v
+            dy_local = -(target.x - cp.x) * sin_v + (target.y - cp.y) * cos_v
+            if 0 < dx_local < self.config.snap_to_target_dist:
+                 offsets.add(dy_local)
+            
+            for offset in offsets:
+                for radius in self.config.sbend_radii:
+                    if abs(offset) >= 2 * radius: continue
+                    self._process_move(current, target, net_width, net_id, open_set, closed_set, snap, f'SB{offset}R{radius}', 'SB', (offset, radius), skip_congestion)
+
+    def _process_move(
+            self,
+            parent: AStarNode,
+            target: Port,
+            net_width: float,
+            net_id: str,
+            open_set: list[AStarNode],
+            closed_set: dict[tuple[int, int, int], float],
+            snap: float,
+            move_type: str,
+            move_class: Literal['S', 'B', 'SB'],
+            params: tuple,
+            skip_congestion: bool,
+            skip_static: bool = False,
+            snap_to_grid: bool = True,
+            ) -> None:
+        cp = parent.port
+        base_ori = round(cp.orientation, 2)
+        state_key = (int(round(cp.x / snap)), int(round(cp.y / snap)), int(round(base_ori / 1.0)))
+        
+        abs_key = (state_key, move_class, params, net_width, self.config.bend_collision_type, snap_to_grid)
+        if abs_key in self._move_cache:
+            res = self._move_cache[abs_key]
+            if move_class == 'B': move_radius = params[0]
+            elif move_class == 'SB': move_radius = params[1]
+            else: move_radius = None
+            self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion)
+            return
+
+        rel_key = (base_ori, move_class, params, net_width, self.config.bend_collision_type, self._self_dilation, snap_to_grid)
+        
+        cache_key = (state_key[0], state_key[1], base_ori, move_type, net_width, snap_to_grid)
+        if cache_key in self._hard_collision_set:
+            return
+            
+        if rel_key in self._move_cache:
+            res_rel = self._move_cache[rel_key]
+            ex = res_rel.end_port.x + cp.x
+            ey = res_rel.end_port.y + cp.y
+            end_state = (int(round(ex / snap)), int(round(ey / snap)), int(round(res_rel.end_port.orientation / 1.0)))
+            if end_state in closed_set and closed_set[end_state] <= parent.g_cost + 1e-6:
+                return
+            res = res_rel.translate(cp.x, cp.y)
+        else:
+            try:
+                if move_class == 'S':
+                    res_rel = Straight.generate(Port(0, 0, base_ori), params[0], net_width, dilation=self._self_dilation, snap_to_grid=snap_to_grid, snap_size=self.config.snap_size)
+                elif move_class == 'B':
+                    res_rel = Bend90.generate(Port(0, 0, base_ori), params[0], net_width, params[1], collision_type=self.config.bend_collision_type, clip_margin=self.config.bend_clip_margin, dilation=self._self_dilation, snap_to_grid=snap_to_grid, snap_size=self.config.snap_size)
+                elif move_class == 'SB':
+                    res_rel = SBend.generate(Port(0, 0, base_ori), params[0], params[1], net_width, collision_type=self.config.bend_collision_type, clip_margin=self.config.bend_clip_margin, dilation=self._self_dilation, snap_to_grid=snap_to_grid, snap_size=self.config.snap_size)
                 else:
-                    rel_key = (base_ori, 'SB', offset, radius, net_width, self.config.bend_collision_type, 0.0)
-                    if rel_key in self._move_cache:
-                        res_rel = self._move_cache[rel_key]
-                        # Check closed set before translating
-                        ex = res_rel.end_port.x + cp.x
-                        ey = res_rel.end_port.y + cp.y
-                        end_state = (round(ex, 3), round(ey, 3), round(res_rel.end_port.orientation, 2))
-                        if end_state in closed_set:
-                            continue
-                        res = res_rel.translate(cp.x, cp.y)
-                    else:
-                        try:
-                            res_rel = SBend.generate(
-                                Port(0, 0, base_ori),
-                                offset,
-                                radius,
-                                width=net_width,
-                                collision_type=self.config.bend_collision_type,
-                                clip_margin=self.config.bend_clip_margin,
-                                dilation=0.0
-                            )
-                            self._move_cache[rel_key] = res_rel
-                            res = res_rel.translate(cp.x, cp.y)
-                        except ValueError:
-                            continue
-                    self._move_cache[abs_key] = res
-                self._add_node(current, res, target, net_width, net_id, open_set, closed_set, f'SB{offset}R{radius}', move_radius=radius)
+                    return
+                self._move_cache[rel_key] = res_rel
+                res = res_rel.translate(cp.x, cp.y)
+            except (ValueError, ZeroDivisionError):
+                return
+
+        self._move_cache[abs_key] = res
+        if move_class == 'B': move_radius = params[0]
+        elif move_class == 'SB': move_radius = params[1]
+        else: move_radius = None
+            
+        self._add_node(parent, res, target, net_width, net_id, open_set, closed_set, move_type, move_radius=move_radius, snap=snap, skip_congestion=skip_congestion)
 
     def _add_node(
             self,
@@ -410,134 +394,83 @@ class AStarRouter:
             net_width: float,
             net_id: str,
             open_set: list[AStarNode],
-            closed_set: set[tuple[float, float, float]],
+            closed_set: dict[tuple[int, int, int], float],
             move_type: str,
             move_radius: float | None = None,
+            snap: float = 1.0,
+            skip_congestion: bool = False,
             ) -> None:
-        # Check closed set before adding to open set
-        state = (round(result.end_port.x, 3), round(result.end_port.y, 3), round(result.end_port.orientation, 2))
-        if state in closed_set:
+        self.metrics['moves_generated'] += 1
+        end_p = result.end_port
+        state = (int(round(end_p.x / snap)), int(round(end_p.y / snap)), int(round(end_p.orientation / 1.0)))
+        
+        if state in closed_set and closed_set[state] <= parent.g_cost + 1e-6:
+            self.metrics['pruned_closed_set'] += 1
             return
 
-        cache_key = (
-            round(parent.port.x, 3),
-            round(parent.port.y, 3),
-            round(parent.port.orientation, 2),
-            move_type,
-            net_width,
-            net_id,
-        )
-        if cache_key in self._collision_cache:
-            if self._collision_cache[cache_key]:
-                return
-        else:
-            # Lazy Dilation: compute dilated polygons only if we need a collision check
-            if result.dilated_geometry is None:
-                # We need to update the ComponentResult with dilated geometry
-                # For simplicity, we'll just buffer the polygons here.
-                # In a more optimized version, ComponentResult might have a .dilate() method.
-                dilated = [p.buffer(self._self_dilation) for p in result.geometry]
-                result.dilated_geometry = dilated
-                # Re-calculate dilated bounds
-                import shapely
-                result.dilated_bounds = shapely.bounds(dilated)
+        parent_p = parent.port
+        cache_key = (int(round(parent_p.x / snap)), int(round(parent_p.y / snap)), int(round(parent_p.orientation / 1.0)), move_type, net_width)
+        
+        if cache_key in self._hard_collision_set:
+            self.metrics['pruned_hard_collision'] += 1
+            return
 
-            hard_coll = False
-            for i, poly in enumerate(result.geometry):
-                dil_poly = result.dilated_geometry[i]
-                if self.cost_evaluator.collision_engine.check_collision(
-                    poly, net_id, buffer_mode='static', start_port=parent.port, end_port=result.end_port,
-                    dilated_geometry=dil_poly
-                ):
-                    hard_coll = True
-                    break
-            self._collision_cache[cache_key] = hard_coll
-            if hard_coll:
-                return
+        is_static_safe = (cache_key in self._static_safe_cache)
+        if not is_static_safe:
+            collision_engine = self.cost_evaluator.collision_engine
+            # Fast check for straights
+            if 'S' in move_type and 'SB' not in move_type:
+                 if collision_engine.check_move_straight_static(parent_p, result.length):
+                      self._hard_collision_set.add(cache_key)
+                      self.metrics['pruned_hard_collision'] += 1
+                      return
+                 is_static_safe = True
+            
+            if not is_static_safe:
+                if collision_engine.check_move_static(result, start_port=parent_p, end_port=end_p):
+                    self._hard_collision_set.add(cache_key)
+                    self.metrics['pruned_hard_collision'] += 1
+                    return
+                else:
+                    self._static_safe_cache.add(cache_key)
 
-        # Lazy Dilation for self-intersection and cost evaluation
-        if result.dilated_geometry is None:
-            dilated = [p.buffer(self._self_dilation) for p in result.geometry]
-            result.dilated_geometry = dilated
-            import shapely
-            result.dilated_bounds = shapely.bounds(dilated)
+        total_overlaps = 0
+        if not skip_congestion:
+            if cache_key in self._congestion_cache:
+                total_overlaps = self._congestion_cache[cache_key]
+            else:
+                total_overlaps = self.cost_evaluator.collision_engine.check_move_congestion(result, net_id)
+                self._congestion_cache[cache_key] = total_overlaps
 
-        # 3. Check for Self-Intersection (Limited to last 100 segments for performance)
-        if result.dilated_geometry:
-            # Union of current move's bounds for fast path-wide pruning
-            m_minx, m_miny, m_maxx, m_maxy = 1e15, 1e15, -1e15, -1e15
-            for b in result.dilated_bounds if result.dilated_bounds is not None else result.bounds:
-
-                m_minx = min(m_minx, b[0])
-                m_miny = min(m_miny, b[1])
-                m_maxx = max(m_maxx, b[2])
-                m_maxy = max(m_maxy, b[3])
-
-            # If current move doesn't overlap the entire parent path bbox, we can skip individual checks
-            # (Except the immediate parent which we usually skip anyway)
-            if parent.path_bbox and not (m_minx > parent.path_bbox[2] or
-                                         m_maxx < parent.path_bbox[0] or
-                                         m_miny > parent.path_bbox[3] or
-                                         m_maxy < parent.path_bbox[1]):
-                
-                for dm_idx, dilated_move in enumerate(result.dilated_geometry):
-                    dm_bounds = result.dilated_bounds[dm_idx]
-                    curr_p: AStarNode | None = parent
-                    seg_idx = 0
-                    while curr_p and curr_p.component_result and seg_idx < 100:
-                        # Skip immediate parent to avoid tangent/port-safety issues
-                        if seg_idx > 0:
-                            res_p = curr_p.component_result
-                            if res_p.dilated_geometry:
-                                for dp_idx, dilated_prev in enumerate(res_p.dilated_geometry):
-                                    dp_bounds = res_p.dilated_bounds[dp_idx]
-                                    # Quick bounds overlap check
-                                    if not (dm_bounds[0] > dp_bounds[2] or
-                                            dm_bounds[2] < dp_bounds[0] or
-                                            dm_bounds[1] > dp_bounds[3] or
-                                            dm_bounds[3] < dp_bounds[1]):
-                                        # Use intersects() which is much faster than intersection()
-                                        if dilated_move.intersects(dilated_prev):
-                                            # Only do expensive area check if absolutely necessary
-                                            overlap = dilated_move.intersection(dilated_prev)
-                                            if not overlap.is_empty and overlap.area > 1e-6:
-                                                return
-                        curr_p = curr_p.parent
-                        seg_idx += 1
+        penalty = 0.0
+        if 'SB' in move_type: penalty = self.config.sbend_penalty
+        elif 'B' in move_type: penalty = self.config.bend_penalty
 
+        # Scale penalty by radius (larger radius = smoother = lower penalty)
+        if move_radius is not None and move_radius > 1e-6:
+            penalty *= (10.0 / move_radius)**0.5
 
         move_cost = self.cost_evaluator.evaluate_move(
-            result.geometry,
-            result.end_port,
-            net_width,
-            net_id,
-            start_port=parent.port,
-            length=result.length,
-            dilated_geometry=result.dilated_geometry,
-            skip_static=True
+            result.geometry, result.end_port, net_width, net_id,
+            start_port=parent_p, length=result.length,
+            dilated_geometry=result.dilated_geometry, penalty=penalty,
+            skip_static=True, skip_congestion=True
         )
+        move_cost += total_overlaps * self.cost_evaluator.congestion_penalty
 
         if move_cost > 1e12:
+            self.metrics['pruned_cost'] += 1
             return
 
-        # Turn penalties scaled by radius to favor larger turns
-        ref_radius = 10.0
-        if 'B' in move_type and move_radius is not None:
-            penalty_factor = ref_radius / move_radius
-            move_cost += self.config.bend_penalty * penalty_factor
-        elif 'SB' in move_type and move_radius is not None:
-            penalty_factor = ref_radius / move_radius
-            move_cost += self.config.sbend_penalty * penalty_factor
-        elif 'B' in move_type:
-            move_cost += self.config.bend_penalty
-        elif 'SB' in move_type:
-            move_cost += self.config.sbend_penalty
-
         g_cost = parent.g_cost + move_cost
+        if state in closed_set and closed_set[state] <= g_cost + 1e-6:
+            self.metrics['pruned_closed_set'] += 1
+            return
+            
         h_cost = self.cost_evaluator.h_manhattan(result.end_port, target)
-
         new_node = AStarNode(result.end_port, g_cost, h_cost, parent, result)
         heapq.heappush(open_set, new_node)
+        self.metrics['moves_added'] += 1
 
     def _reconstruct_path(self, end_node: AStarNode) -> list[ComponentResult]:
         path = []

inire/router/config.py

@@ -10,13 +10,24 @@ class RouterConfig:
     """Configuration parameters for the A* Router."""
 
     node_limit: int = 1000000
-    straight_lengths: list[float] = field(default_factory=lambda: [1.0, 5.0, 25.0])
-    bend_radii: list[float] = field(default_factory=lambda: [10.0])
-    sbend_offsets: list[float] = field(default_factory=lambda: [-5.0, -2.0, 2.0, 5.0])
-    sbend_radii: list[float] = field(default_factory=lambda: [10.0])
-    snap_to_target_dist: float = 20.0
-    bend_penalty: float = 50.0
-    sbend_penalty: float = 100.0
+    snap_size: float = 5.0
+    # Sparse Sampling Configuration
+    max_straight_length: float = 2000.0
+    num_straight_samples: int = 5
+    min_straight_length: float = 5.0
+    
+    # Offsets for SBends (still list-based for now, or could range)
+    sbend_offsets: list[float] = field(default_factory=lambda: [-100.0, -50.0, -10.0, 10.0, 50.0, 100.0])
+    
+    # Deprecated but kept for compatibility during refactor
+    straight_lengths: list[float] = field(default_factory=list) 
+    
+    bend_radii: list[float] = field(default_factory=lambda: [50.0, 100.0])
+    sbend_radii: list[float] = field(default_factory=lambda: [50.0, 100.0, 500.0])
+    snap_to_target_dist: float = 1000.0
+    use_analytical_sbends: bool = True
+    bend_penalty: float = 250.0
+    sbend_penalty: float = 500.0
     bend_collision_type: Literal["arc", "bbox", "clipped_bbox"] | Any = "arc"
     bend_clip_margin: float = 10.0
 
@@ -28,4 +39,5 @@ class CostConfig:
     unit_length_cost: float = 1.0
     greedy_h_weight: float = 1.1
     congestion_penalty: float = 10000.0
-    bend_penalty: float = 50.0
+    bend_penalty: float = 250.0
+    sbend_penalty: float = 500.0

inire/router/cost.py

@@ -39,7 +39,8 @@ class CostEvaluator:
             unit_length_cost: float = 1.0,
             greedy_h_weight: float = 1.1,
             congestion_penalty: float = 10000.0,
-            bend_penalty: float = 50.0,
+            bend_penalty: float = 250.0,
+            sbend_penalty: float = 500.0,
             ) -> None:
         """
         Initialize the Cost Evaluator.
@@ -51,6 +52,7 @@ class CostEvaluator:
             greedy_h_weight: Heuristic weighting (A* greedy factor).
             congestion_penalty: Multiplier for path overlaps in negotiated congestion.
             bend_penalty: Base cost for 90-degree bends.
+            sbend_penalty: Base cost for parametric S-bends.
         """
         self.collision_engine = collision_engine
         self.danger_map = danger_map
@@ -59,6 +61,7 @@ class CostEvaluator:
             greedy_h_weight=greedy_h_weight,
             congestion_penalty=congestion_penalty,
             bend_penalty=bend_penalty,
+            sbend_penalty=sbend_penalty,
         )
 
         # Use config values
@@ -81,29 +84,18 @@ class CostEvaluator:
 
     def h_manhattan(self, current: Port, target: Port) -> float:
         """
-        Heuristic: weighted Manhattan distance + orientation penalty.
-
-        Args:
-            current: Current port state.
-            target: Target port state.
-
-        Returns:
-            Heuristic cost estimate.
+        Heuristic: weighted Manhattan distance + mandatory turn penalties.
         """
         dx = abs(current.x - target.x)
         dy = abs(current.y - target.y)
         dist = dx + dy
 
-        # Orientation penalty if not aligned with target entry
-        # If we need to turn, the cost is at least min_bend_radius * pi/2
-        # But we also need to account for the physical distance required for the turn.
         penalty = 0.0
-        if current.orientation != target.orientation:
-            # 90-degree turn cost: radius 10 -> ~15.7 um + penalty
-            penalty += 15.7 + self.config.bend_penalty
+        if abs(current.orientation - target.orientation) > 0.1:
+            # Needs at least 1 bend
+            penalty += 10.0 + self.config.bend_penalty * 0.1
 
-        # Add 1.5 multiplier for greediness (faster search)
-        return 1.5 * (dist + penalty)
+        return self.greedy_h_weight * (dist + penalty)
 
 
     def evaluate_move(
@@ -116,6 +108,8 @@ class CostEvaluator:
             length: float = 0.0,
             dilated_geometry: list[Polygon] | None = None,
             skip_static: bool = False,
+            skip_congestion: bool = False,
+            penalty: float = 0.0,
             ) -> float:
         """
         Calculate the cost of a single move (Straight, Bend, SBend).
@@ -128,36 +122,44 @@ class CostEvaluator:
             start_port: Port at the start of the move.
             length: Physical path length of the move.
             dilated_geometry: Pre-calculated dilated polygons.
-            skip_static: If True, bypass static collision checks (e.g. if already done).
+            skip_static: If True, bypass static collision checks.
+            skip_congestion: If True, bypass congestion checks.
+            penalty: Fixed cost penalty for the move type.
 
         Returns:
             Total cost of the move, or 1e15 if invalid.
         """
         _ = net_width  # Unused
-        total_cost = length * self.unit_length_cost
-
+        
         # 1. Boundary Check
-        if not self.danger_map.is_within_bounds(end_port.x, end_port.y):
+        danger_map = self.danger_map
+        if not danger_map.is_within_bounds(end_port.x, end_port.y):
             return 1e15
 
-        # 2. Collision Check
-        for i, poly in enumerate(geometry):
-            dil_poly = dilated_geometry[i] if dilated_geometry else None
-            # Hard Collision (Static obstacles)
-            if not skip_static:
-                if self.collision_engine.check_collision(
-                    poly, net_id, buffer_mode='static', start_port=start_port, end_port=end_port,
-                    dilated_geometry=dil_poly
-                ):
-                    return 1e15
+        total_cost = length * self.unit_length_cost + penalty
 
-            # Soft Collision (Negotiated Congestion)
-            overlaps = self.collision_engine.check_collision(
-                poly, net_id, buffer_mode='congestion', dilated_geometry=dil_poly
-            )
-            if isinstance(overlaps, int) and overlaps > 0:
-                total_cost += overlaps * self.congestion_penalty
+        # 2. Collision Check
+        # FAST PATH: skip_static and skip_congestion are often True when called from optimized AStar
+        if not skip_static or not skip_congestion:
+            collision_engine = self.collision_engine
+            for i, poly in enumerate(geometry):
+                dil_poly = dilated_geometry[i] if dilated_geometry else None
+                # Hard Collision (Static obstacles)
+                if not skip_static:
+                    if collision_engine.check_collision(
+                        poly, net_id, buffer_mode='static', start_port=start_port, end_port=end_port,
+                        dilated_geometry=dil_poly
+                    ):
+                        return 1e15
+
+                # Soft Collision (Negotiated Congestion)
+                if not skip_congestion:
+                    overlaps = collision_engine.check_collision(
+                        poly, net_id, buffer_mode='congestion', dilated_geometry=dil_poly
+                    )
+                    if isinstance(overlaps, int) and overlaps > 0:
+                        total_cost += overlaps * self.congestion_penalty
 
         # 3. Proximity cost from Danger Map
-        total_cost += self.g_proximity(end_port.x, end_port.y)
+        total_cost += danger_map.get_cost(end_port.x, end_port.y)
         return total_cost

inire/router/pathfinder.py

@@ -2,8 +2,9 @@ from __future__ import annotations
 
 import logging
 import time
+import random
 from dataclasses import dataclass
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Callable
 
 if TYPE_CHECKING:
     from inire.geometry.components import ComponentResult
@@ -26,17 +27,20 @@ class RoutingResult:
     """ List of moves forming the path """
 
     is_valid: bool
-    """ Whether the path is collision-free """
+    """ Whether the path is collision-free and reached the target """
 
     collisions: int
     """ Number of detected collisions/overlaps """
 
+    reached_target: bool = False
+    """ Whether the final port matches the target port """
+
 
 class PathFinder:
     """
     Multi-net router using Negotiated Congestion.
     """
-    __slots__ = ('router', 'cost_evaluator', 'max_iterations', 'base_congestion_penalty')
+    __slots__ = ('router', 'cost_evaluator', 'max_iterations', 'base_congestion_penalty', 'use_tiered_strategy', 'congestion_multiplier', 'accumulated_expanded_nodes')
 
     router: AStarRouter
     """ The A* search engine """
@@ -50,12 +54,20 @@ class PathFinder:
     base_congestion_penalty: float
     """ Starting penalty for overlaps """
 
+    congestion_multiplier: float
+    """ Multiplier for congestion penalty per iteration """
+
+    use_tiered_strategy: bool
+    """ If True, use simpler collision models in early iterations for speed """
+
     def __init__(
             self,
             router: AStarRouter,
             cost_evaluator: CostEvaluator,
             max_iterations: int = 10,
             base_congestion_penalty: float = 100.0,
+            congestion_multiplier: float = 1.5,
+            use_tiered_strategy: bool = True,
             ) -> None:
         """
         Initialize the PathFinder.
@@ -65,16 +77,25 @@ class PathFinder:
             cost_evaluator: The evaluator for path costs.
             max_iterations: Maximum number of rip-up and reroute iterations.
             base_congestion_penalty: Starting penalty for overlaps.
+            congestion_multiplier: Multiplier for congestion penalty per iteration.
+            use_tiered_strategy: Whether to use simplified collision models in early iterations.
         """
         self.router = router
         self.cost_evaluator = cost_evaluator
         self.max_iterations = max_iterations
         self.base_congestion_penalty = base_congestion_penalty
+        self.congestion_multiplier = congestion_multiplier
+        self.use_tiered_strategy = use_tiered_strategy
+        self.accumulated_expanded_nodes: list[tuple[float, float, float]] = []
 
     def route_all(
             self,
             netlist: dict[str, tuple[Port, Port]],
             net_widths: dict[str, float],
+            store_expanded: bool = False,
+            iteration_callback: Callable[[int, dict[str, RoutingResult]], None] | None = None,
+            shuffle_nets: bool = False,
+            seed: int | None = None,
             ) -> dict[str, RoutingResult]:
         """
         Route all nets in the netlist using Negotiated Congestion.
@@ -82,23 +103,40 @@ class PathFinder:
         Args:
             netlist: Mapping of net_id to (start_port, target_port).
             net_widths: Mapping of net_id to waveguide width.
+            store_expanded: Whether to store expanded nodes for ALL iterations and nets.
+            iteration_callback: Optional callback(iteration_idx, current_results).
+            shuffle_nets: Whether to randomize the order of nets each iteration.
+            seed: Optional seed for randomization (enables reproducibility).
 
         Returns:
             Mapping of net_id to RoutingResult.
         """
         results: dict[str, RoutingResult] = {}
         self.cost_evaluator.congestion_penalty = self.base_congestion_penalty
+        self.accumulated_expanded_nodes = []
 
         start_time = time.monotonic()
         num_nets = len(netlist)
         session_timeout = max(60.0, 10.0 * num_nets * self.max_iterations)
 
+        all_net_ids = list(netlist.keys())
+
         for iteration in range(self.max_iterations):
             any_congestion = False
+            # Clear accumulation for this iteration so callback gets fresh data
+            self.accumulated_expanded_nodes = []
+            
             logger.info(f'PathFinder Iteration {iteration}...')
 
+            # 0. Shuffle nets if requested
+            if shuffle_nets:
+                # Use a new seed based on iteration for deterministic different orders
+                it_seed = (seed + iteration) if seed is not None else None
+                random.Random(it_seed).shuffle(all_net_ids)
+
             # Sequence through nets
-            for net_id, (start, target) in netlist.items():
+            for net_id in all_net_ids:
+                start, target = netlist[net_id]
                 # Timeout check
                 elapsed = time.monotonic() - start_time
                 if elapsed > session_timeout:
@@ -111,51 +149,128 @@ class PathFinder:
                 self.cost_evaluator.collision_engine.remove_path(net_id)
 
                 # 2. Reroute with current congestion info
-                # Tiered Strategy: use clipped_bbox for Iteration 0 for speed.
-                # Switch to arc for higher iterations if collisions persist.
-                coll_model = "clipped_bbox" if iteration == 0 else "arc"
+                target_coll_model = self.router.config.bend_collision_type
+                coll_model = target_coll_model
+                skip_cong = False
+                if self.use_tiered_strategy and iteration == 0:
+                    skip_cong = True
+                    if target_coll_model == "arc":
+                        coll_model = "clipped_bbox"
+                
+                # Dynamic node limit: increase if it failed previously
+                base_node_limit = self.router.config.node_limit
+                current_node_limit = base_node_limit
+                if net_id in results and not results[net_id].reached_target:
+                     current_node_limit = base_node_limit * (iteration + 1)
                 
                 net_start = time.monotonic()
-                path = self.router.route(start, target, width, net_id=net_id, bend_collision_type=coll_model)
+                
+                # Temporarily override node_limit
+                original_limit = self.router.node_limit
+                self.router.node_limit = current_node_limit
+                
+                path = self.router.route(start, target, width, net_id=net_id, bend_collision_type=coll_model, return_partial=True, store_expanded=store_expanded, skip_congestion=skip_cong)
+                
+                if store_expanded and self.router.last_expanded_nodes:
+                    self.accumulated_expanded_nodes.extend(self.router.last_expanded_nodes)
+
+                # Restore
+                self.router.node_limit = original_limit
+                
                 logger.debug(f'  Net {net_id} routed in {time.monotonic() - net_start:.4f}s using {coll_model}')
 
                 if path:
-                    # 3. Add to index
+                    # Check if reached exactly
+                    last_p = path[-1].end_port
+                    reached = (abs(last_p.x - target.x) < 1e-6 and
+                               abs(last_p.y - target.y) < 1e-6 and
+                               abs(last_p.orientation - target.orientation) < 0.1)
+                    
                     all_geoms = []
                     all_dilated = []
-                    for res in path:
-                        all_geoms.extend(res.geometry)
-                        if res.dilated_geometry:
-                            all_dilated.extend(res.dilated_geometry)
-                        else:
-                            # Fallback dilation
-                            dilation = self.cost_evaluator.collision_engine.clearance / 2.0
-                            all_dilated.extend([p.buffer(dilation) for p in res.geometry])
+                    # 3. Add to index ONLY if it reached the target
+                    # (Prevents failed paths from blocking others forever)
+                    if reached:
+                        for res in path:
+                            # Use the search geometry (could be proxy or arc) for indexing
+                            # to ensure consistency with what other nets use for their search.
+                            all_geoms.extend(res.geometry)
                             
-                    self.cost_evaluator.collision_engine.add_path(net_id, all_geoms, dilated_geometry=all_dilated)
+                            if res.dilated_geometry:
+                                all_dilated.extend(res.dilated_geometry)
+                            else:
+                                # Fallback dilation
+                                dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+                                all_dilated.extend([p.buffer(dilation) for p in res.geometry])
+                                
+                        self.cost_evaluator.collision_engine.add_path(net_id, all_geoms, dilated_geometry=all_dilated)
 
                     # Check if this new path has any congestion
                     collision_count = 0
-                    for i, poly in enumerate(all_geoms):
-                        overlaps = self.cost_evaluator.collision_engine.check_collision(
-                            poly, net_id, buffer_mode='congestion', dilated_geometry=all_dilated[i]
-                        )
-                        if isinstance(overlaps, int):
-                            collision_count += overlaps
+                    # Always check for congestion to decide if more iterations are needed
+                    if reached:
+                        # For FINAL verification of this net's success, we should ideally
+                        # use high-fidelity geometry if available, but since Negotiated
+                        # Congestion relies on what is IN the index, we check the indexed geoms.
+                        # BUT, to fix the "false failed" issue where clipped_bbox overlaps
+                        # even if arcs don't, we should verify with actual_geometry.
+                        
+                        verif_geoms = []
+                        verif_dilated = []
+                        for res in path:
+                            is_proxy = (res.actual_geometry is not None)
+                            g = res.actual_geometry if is_proxy else res.geometry
+                            verif_geoms.extend(g)
+                            
+                            # If we are using actual_geometry as high-fidelity replacement for a proxy,
+                            # we MUST ensure we use the high-fidelity dilation too.
+                            if is_proxy:
+                                # ComponentResult stores dilated_geometry for the 'geometry' (proxy).
+                                # It does NOT store it for 'actual_geometry' unless we re-buffer.
+                                dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+                                verif_dilated.extend([p.buffer(dilation) for p in g])
+                            else:
+                                # Use existing dilated geometry if it matches the current geom
+                                if res.dilated_geometry:
+                                    verif_dilated.extend(res.dilated_geometry)
+                                else:
+                                    dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+                                    verif_dilated.extend([p.buffer(dilation) for p in g])
+
+                        for i, poly in enumerate(verif_geoms):
+                            # IMPORTANT: We check against OTHER nets. 
+                            # If we just check self.check_congestion(poly, net_id), 
+                            # it checks against the dynamic index which ALREADY contains this net's 
+                            # path (added in step 3 above).
+                            # To correctly count REAL overlaps with others:
+                            self.cost_evaluator.collision_engine._ensure_dynamic_tree()
+                            if self.cost_evaluator.collision_engine.dynamic_tree:
+                                hits = self.cost_evaluator.collision_engine.dynamic_tree.query(verif_dilated[i], predicate='intersects')
+                                for hit_idx in hits:
+                                    obj_id = self.cost_evaluator.collision_engine.dynamic_obj_ids[hit_idx]
+                                    other_net_id, _ = self.cost_evaluator.collision_engine.dynamic_geometries[obj_id]
+                                    if other_net_id != net_id:
+                                        collision_count += 1
 
                     if collision_count > 0:
                         any_congestion = True
 
-                    results[net_id] = RoutingResult(net_id, path, collision_count == 0, collision_count)
+                    results[net_id] = RoutingResult(net_id, path, (reached and collision_count == 0), collision_count, reached_target=reached)
                 else:
-                    results[net_id] = RoutingResult(net_id, [], False, 0)
+                    results[net_id] = RoutingResult(net_id, [], False, 0, reached_target=False)
                     any_congestion = True
 
+            if iteration_callback:
+                iteration_callback(iteration, results)
+
             if not any_congestion:
-                break
+                # Check if all reached target
+                all_reached = all(r.reached_target for r in results.values())
+                if all_reached:
+                    break
 
             # 4. Inflate congestion penalty
-            self.cost_evaluator.congestion_penalty *= 1.5
+            self.cost_evaluator.congestion_penalty *= self.congestion_multiplier
 
         return self._finalize_results(results, netlist)
 
@@ -183,15 +298,41 @@ class PathFinder:
                 continue
 
             collision_count = 0
+            # Use high-fidelity verification against OTHER nets
+            verif_geoms = []
+            verif_dilated = []
             for comp in res.path:
-                for i, poly in enumerate(comp.geometry):
-                    dil_poly = comp.dilated_geometry[i] if comp.dilated_geometry else None
-                    overlaps = self.cost_evaluator.collision_engine.check_collision(
-                        poly, net_id, buffer_mode='congestion', dilated_geometry=dil_poly
-                    )
-                    if isinstance(overlaps, int):
-                        collision_count += overlaps
+                is_proxy = (comp.actual_geometry is not None)
+                g = comp.actual_geometry if is_proxy else comp.geometry
+                verif_geoms.extend(g)
+                if is_proxy:
+                    dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+                    verif_dilated.extend([p.buffer(dilation) for p in g])
+                else:
+                    if comp.dilated_geometry:
+                        verif_dilated.extend(comp.dilated_geometry)
+                    else:
+                        dilation = self.cost_evaluator.collision_engine.clearance / 2.0
+                        verif_dilated.extend([p.buffer(dilation) for p in g])
 
-            final_results[net_id] = RoutingResult(net_id, res.path, collision_count == 0, collision_count)
+            self.cost_evaluator.collision_engine._ensure_dynamic_tree()
+            if self.cost_evaluator.collision_engine.dynamic_tree:
+                for i, poly in enumerate(verif_geoms):
+                    hits = self.cost_evaluator.collision_engine.dynamic_tree.query(verif_dilated[i], predicate='intersects')
+                    for hit_idx in hits:
+                        obj_id = self.cost_evaluator.collision_engine.dynamic_obj_ids[hit_idx]
+                        other_net_id, _ = self.cost_evaluator.collision_engine.dynamic_geometries[obj_id]
+                        if other_net_id != net_id:
+                            collision_count += 1
+
+            reached = False
+            if res.path:
+                target_p = netlist[net_id][1]
+                last_p = res.path[-1].end_port
+                reached = (abs(last_p.x - target_p.x) < 1e-6 and
+                           abs(last_p.y - target_p.y) < 1e-6 and
+                           abs(last_p.orientation - target_p.orientation) < 0.1)
+
+            final_results[net_id] = RoutingResult(net_id, res.path, (collision_count == 0 and reached), collision_count, reached_target=reached)
 
         return final_results

inire/router/visibility.py

@@ -0,0 +1,125 @@
+from __future__ import annotations
+
+import numpy
+from typing import TYPE_CHECKING
+import rtree
+from shapely.geometry import Point, LineString
+
+if TYPE_CHECKING:
+    from inire.geometry.collision import CollisionEngine
+    from inire.geometry.primitives import Port
+
+
+from inire.geometry.primitives import Port
+
+class VisibilityManager:
+    """
+    Manages corners of static obstacles for sparse A* / Visibility Graph jumps.
+    """
+    __slots__ = ('collision_engine', 'corners', 'corner_index', '_corner_graph', '_static_visibility_cache')
+
+    def __init__(self, collision_engine: CollisionEngine) -> None:
+        self.collision_engine = collision_engine
+        self.corners: list[tuple[float, float]] = []
+        self.corner_index = rtree.index.Index()
+        self._corner_graph: dict[int, list[tuple[float, float, float]]] = {}
+        self._static_visibility_cache: dict[tuple[int, int], list[tuple[float, float, float]]] = {}
+        self._build()
+
+    def _build(self) -> None:
+        """
+        Extract corners and pre-compute corner-to-corner visibility.
+        """
+        raw_corners = []
+        for obj_id, poly in self.collision_engine.static_dilated.items():
+            coords = list(poly.exterior.coords)
+            if coords[0] == coords[-1]:
+                coords = coords[:-1]
+            raw_corners.extend(coords)
+            for ring in poly.interiors:
+                coords = list(ring.coords)
+                if coords[0] == coords[-1]:
+                    coords = coords[:-1]
+                raw_corners.extend(coords)
+
+        if not raw_corners:
+            return
+
+        # Deduplicate and snap to 1nm
+        seen = set()
+        for x, y in raw_corners:
+            sx, sy = round(x, 3), round(y, 3)
+            if (sx, sy) not in seen:
+                seen.add((sx, sy))
+                self.corners.append((sx, sy))
+
+        # Build spatial index for corners
+        for i, (x, y) in enumerate(self.corners):
+            self.corner_index.insert(i, (x, y, x, y))
+
+        # Pre-compute visibility graph between corners
+        num_corners = len(self.corners)
+        if num_corners > 200:
+             # Limit pre-computation if too many corners
+             return
+
+        for i in range(num_corners):
+            self._corner_graph[i] = []
+            p1 = Port(self.corners[i][0], self.corners[i][1], 0)
+            for j in range(num_corners):
+                if i == j: continue
+                cx, cy = self.corners[j]
+                dx, dy = cx - p1.x, cy - p1.y
+                dist = numpy.sqrt(dx**2 + dy**2)
+                angle = numpy.degrees(numpy.arctan2(dy, dx))
+                reach = self.collision_engine.ray_cast(p1, angle, max_dist=dist + 0.05)
+                if reach >= dist - 0.01:
+                    self._corner_graph[i].append((cx, cy, dist))
+
+    def get_visible_corners(self, origin: Port, max_dist: float = 1000.0) -> list[tuple[float, float, float]]:
+        """
+        Find all corners visible from the origin.
+        Returns list of (x, y, distance).
+        """
+        if max_dist < 0:
+            return []
+            
+        ox, oy = round(origin.x, 3), round(origin.y, 3)
+        
+        # 1. Exact corner check
+        # Use spatial index to find if origin is AT a corner
+        nearby = list(self.corner_index.intersection((ox - 0.001, oy - 0.001, ox + 0.001, oy + 0.001)))
+        for idx in nearby:
+            cx, cy = self.corners[idx]
+            if abs(cx - ox) < 1e-4 and abs(cy - oy) < 1e-4:
+                # We are at a corner! Return pre-computed graph (filtered by max_dist)
+                if idx in self._corner_graph:
+                    return [c for c in self._corner_graph[idx] if c[2] <= max_dist]
+
+        # 2. Cache check for arbitrary points
+        # Grid-based caching for arbitrary points is tricky, 
+        # but since static obstacles don't change, we can cache exact coordinates.
+        cache_key = (int(ox * 1000), int(oy * 1000))
+        if cache_key in self._static_visibility_cache:
+             return self._static_visibility_cache[cache_key]
+
+        # 3. Full visibility check
+        bounds = (origin.x - max_dist, origin.y - max_dist, origin.x + max_dist, origin.y + max_dist)
+        candidates = list(self.corner_index.intersection(bounds))
+        
+        visible = []
+        for i in candidates:
+            cx, cy = self.corners[i]
+            dx, dy = cx - origin.x, cy - origin.y
+            dist = numpy.sqrt(dx**2 + dy**2)
+            
+            if dist > max_dist or dist < 1e-3:
+                continue
+                
+            angle = numpy.degrees(numpy.arctan2(dy, dx))
+            reach = self.collision_engine.ray_cast(origin, angle, max_dist=dist + 0.05)
+            if reach >= dist - 0.01:
+                visible.append((cx, cy, dist))
+        
+        self._static_visibility_cache[cache_key] = visible
+        return visible

inire/tests/test_astar.py

@@ -15,11 +15,11 @@ def basic_evaluator() -> CostEvaluator:
     engine = CollisionEngine(clearance=2.0)
     danger_map = DangerMap(bounds=(0, 0, 100, 100))
     danger_map.precompute([])
-    return CostEvaluator(engine, danger_map)
+    return CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
 
 
 def test_astar_straight(basic_evaluator: CostEvaluator) -> None:
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0])
     start = Port(0, 0, 0)
     target = Port(50, 0, 0)
     path = router.route(start, target, net_width=2.0)
@@ -35,11 +35,9 @@ def test_astar_straight(basic_evaluator: CostEvaluator) -> None:
 
 
 def test_astar_bend(basic_evaluator: CostEvaluator) -> None:
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     start = Port(0, 0, 0)
     # 20um right, 20um up. Needs a 10um bend and a 10um bend.
-    # From (0,0,0) -> Bend90 CW R=10 -> (10, -10, 270) ??? No.
-    # Try: (0,0,0) -> Bend90 CCW R=10 -> (10, 10, 90) -> Straight 10 -> (10, 20, 90) -> Bend90 CW R=10 -> (20, 30, 0)
     target = Port(20, 20, 0)
     path = router.route(start, target, net_width=2.0)
 
@@ -58,7 +56,7 @@ def test_astar_obstacle(basic_evaluator: CostEvaluator) -> None:
     basic_evaluator.collision_engine.add_static_obstacle(obstacle)
     basic_evaluator.danger_map.precompute([obstacle])
 
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0, straight_lengths=[1.0, 5.0, 25.0], bend_radii=[10.0])
     router.node_limit = 1000000 # Give it more room for detour
     start = Port(0, 0, 0)
     target = Port(60, 0, 0)
@@ -74,7 +72,7 @@ def test_astar_obstacle(basic_evaluator: CostEvaluator) -> None:
 
 
 def test_astar_snap_to_target_lookahead(basic_evaluator: CostEvaluator) -> None:
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0)
     # Target is NOT on 1um grid
     start = Port(0, 0, 0)
     target = Port(10.1, 0, 0)

inire/tests/test_components.py

@@ -8,7 +8,7 @@ def test_straight_generation() -> None:
     start = Port(0, 0, 0)
     length = 10.0
     width = 2.0
-    result = Straight.generate(start, length, width)
+    result = Straight.generate(start, length, width, snap_size=1.0)
 
     assert result.end_port.x == 10.0
     assert result.end_port.y == 0.0
@@ -29,13 +29,13 @@ def test_bend90_generation() -> None:
     width = 2.0
 
     # CW bend
-    result_cw = Bend90.generate(start, radius, width, direction="CW")
+    result_cw = Bend90.generate(start, radius, width, direction="CW", snap_size=1.0)
     assert result_cw.end_port.x == 10.0
     assert result_cw.end_port.y == -10.0
     assert result_cw.end_port.orientation == 270.0
 
     # CCW bend
-    result_ccw = Bend90.generate(start, radius, width, direction="CCW")
+    result_ccw = Bend90.generate(start, radius, width, direction="CCW", snap_size=1.0)
     assert result_ccw.end_port.x == 10.0
     assert result_ccw.end_port.y == 10.0
     assert result_ccw.end_port.orientation == 90.0
@@ -47,7 +47,7 @@ def test_sbend_generation() -> None:
     radius = 10.0
     width = 2.0
 
-    result = SBend.generate(start, offset, radius, width)
+    result = SBend.generate(start, offset, radius, width, snap_size=1.0)
     assert result.end_port.y == 5.0
     assert result.end_port.orientation == 0.0
     assert len(result.geometry) == 2 # Optimization: returns individual arcs
@@ -63,7 +63,7 @@ def test_bend_collision_models() -> None:
     width = 2.0
 
     # 1. BBox model
-    res_bbox = Bend90.generate(start, radius, width, direction="CCW", collision_type="bbox")
+    res_bbox = Bend90.generate(start, radius, width, direction="CCW", collision_type="bbox", snap_size=1.0)
     # Arc CCW R=10 from (0,0,0) ends at (10,10,90).
     # Waveguide width is 2.0, so bbox will be slightly larger than (0,0,10,10)
     minx, miny, maxx, maxy = res_bbox.geometry[0].bounds
@@ -73,7 +73,7 @@ def test_bend_collision_models() -> None:
     assert maxy >= 10.0 - 1e-6
 
     # 2. Clipped BBox model
-    res_clipped = Bend90.generate(start, radius, width, direction="CCW", collision_type="clipped_bbox", clip_margin=1.0)
+    res_clipped = Bend90.generate(start, radius, width, direction="CCW", collision_type="clipped_bbox", clip_margin=1.0, snap_size=1.0)
     # Area should be less than full bbox
     assert res_clipped.geometry[0].area < res_bbox.geometry[0].area
 
@@ -84,11 +84,11 @@ def test_sbend_collision_models() -> None:
     radius = 10.0
     width = 2.0
 
-    res_bbox = SBend.generate(start, offset, radius, width, collision_type="bbox")
+    res_bbox = SBend.generate(start, offset, radius, width, collision_type="bbox", snap_size=1.0)
     # Geometry should be a list of individual bbox polygons for each arc
     assert len(res_bbox.geometry) == 2
 
-    res_arc = SBend.generate(start, offset, radius, width, collision_type="arc")
+    res_arc = SBend.generate(start, offset, radius, width, collision_type="arc", snap_size=1.0)
     area_bbox = sum(p.area for p in res_bbox.geometry)
     area_arc = sum(p.area for p in res_arc.geometry)
     assert area_bbox > area_arc
@@ -101,7 +101,8 @@ def test_sbend_continuity() -> None:
     radius = 20.0
     width = 1.0
 
-    res = SBend.generate(start, offset, radius, width)
+    # We use snap_size=1.0 so that (10-offset) = 6.0 is EXACTLY hit.
+    res = SBend.generate(start, offset, radius, width, snap_size=1.0)
 
     # Target orientation should be same as start
     assert abs(res.end_port.orientation - 90.0) < 1e-6
@@ -122,9 +123,9 @@ def test_arc_sagitta_precision() -> None:
     width = 2.0
 
     # Coarse: 1um sagitta
-    res_coarse = Bend90.generate(start, radius, width, sagitta=1.0)
+    res_coarse = Bend90.generate(start, radius, width, direction="CCW", sagitta=1.0)
     # Fine: 0.01um (10nm) sagitta
-    res_fine = Bend90.generate(start, radius, width, sagitta=0.01)
+    res_fine = Bend90.generate(start, radius, width, direction="CCW", sagitta=0.01)
 
     # Number of segments should be significantly higher for fine
     # Exterior points = (segments + 1) * 2
@@ -141,7 +142,7 @@ def test_component_transform_invariance() -> None:
     radius = 10.0
     width = 2.0
 
-    res0 = Bend90.generate(start0, radius, width, direction="CCW")
+    res0 = Bend90.generate(start0, radius, width, direction="CCW", snap_size=1.0)
 
     # Transform: Translate (10, 10) then Rotate 90
     dx, dy = 10.0, 5.0
@@ -152,7 +153,7 @@ def test_component_transform_invariance() -> None:
 
     # 2. Generate at transformed start
     start_transformed = rotate_port(translate_port(start0, dx, dy), angle)
-    res_transformed = Bend90.generate(start_transformed, radius, width, direction="CCW")
+    res_transformed = Bend90.generate(start_transformed, radius, width, direction="CCW", snap_size=1.0)
 
     assert abs(res_transformed.end_port.x - p_end_transformed.x) < 1e-6
     assert abs(res_transformed.end_port.y - p_end_transformed.y) < 1e-6

inire/tests/test_congestion.py

@@ -15,11 +15,11 @@ def basic_evaluator() -> CostEvaluator:
     # Wider bounds to allow going around (y from -40 to 40)
     danger_map = DangerMap(bounds=(0, -40, 100, 40))
     danger_map.precompute([])
-    return CostEvaluator(engine, danger_map)
+    return CostEvaluator(engine, danger_map, bend_penalty=50.0, sbend_penalty=150.0)
 
 
 def test_astar_sbend(basic_evaluator: CostEvaluator) -> None:
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0, sbend_offsets=[2.0, 5.0])
     # Start at (0,0), target at (50, 2) -> 2um lateral offset
     # This matches one of our discretized SBend offsets.
     start = Port(0, 0, 0)
@@ -39,7 +39,7 @@ def test_astar_sbend(basic_evaluator: CostEvaluator) -> None:
 
 
 def test_pathfinder_negotiated_congestion_resolution(basic_evaluator: CostEvaluator) -> None:
-    router = AStarRouter(basic_evaluator)
+    router = AStarRouter(basic_evaluator, snap_size=1.0, bend_radii=[5.0, 10.0], sbend_radii=[5.0, 10.0])
     # Increase base penalty to force detour immediately
     pf = PathFinder(router, basic_evaluator, max_iterations=10, base_congestion_penalty=1000.0)
 

inire/tests/test_refinements.py

@@ -11,9 +11,10 @@ def test_arc_resolution_sagitta() -> None:
     start = Port(0, 0, 0)
     # R=10, 90 deg bend.
     # High tolerance (0.5um) -> few segments
-    res_coarse = Bend90.generate(start, radius=10.0, width=2.0, sagitta=0.5)
-    # Low tolerance (0.001um = 1nm) -> many segments
-    res_fine = Bend90.generate(start, radius=10.0, width=2.0, sagitta=0.001)
+    res_coarse = Bend90.generate(start, radius=10.0, width=2.0, direction="CCW", sagitta=0.5)
+    # Low tolerance (1nm) -> many segments
+    res_fine = Bend90.generate(start, radius=10.0, width=2.0, direction="CCW", sagitta=0.001)
+
 
     # Check number of points in the polygon exterior
     # (num_segments + 1) * 2 points usually
@@ -28,7 +29,7 @@ def test_locked_paths() -> None:
     danger_map = DangerMap(bounds=(0, -50, 100, 50))
     danger_map.precompute([])
     evaluator = CostEvaluator(engine, danger_map)
-    router = AStarRouter(evaluator)
+    router = AStarRouter(evaluator, bend_radii=[5.0, 10.0], sbend_radii=[5.0, 10.0])
     pf = PathFinder(router, evaluator)
 
     # 1. Route Net A

inire/utils/visualization.py

@@ -18,6 +18,7 @@ def plot_routing_results(
     static_obstacles: list[Polygon],
     bounds: tuple[float, float, float, float],
     netlist: dict[str, tuple[Port, Port]] | None = None,
+    show_actual: bool = True,
 ) -> tuple[Figure, Axes]:
     """
     Plot obstacles and routed paths using matplotlib.
@@ -27,62 +28,191 @@ def plot_routing_results(
         static_obstacles: List of static obstacle polygons.
         bounds: Plot limits (minx, miny, maxx, maxy).
         netlist: Optional original netlist for port visualization.
+        show_actual: If True, overlay high-fidelity geometry if available.
 
     Returns:
         The matplotlib Figure and Axes objects.
     """
-    fig, ax = plt.subplots(figsize=(10, 10))
+    fig, ax = plt.subplots(figsize=(12, 12))
 
     # Plot static obstacles (gray)
     for poly in static_obstacles:
         x, y = poly.exterior.xy
-        ax.fill(x, y, alpha=0.5, fc="gray", ec="black")
+        ax.fill(x, y, alpha=0.3, fc="gray", ec="black", zorder=1)
 
     # Plot paths
-    colors = plt.get_cmap("tab10")
+    colors = plt.get_cmap("tab20")
     for i, (net_id, res) in enumerate(results.items()):
-        # Use modulo to avoid index out of range for many nets
-        color: str | tuple[float, ...] = colors(i % 10)
+        color: str | tuple[float, ...] = colors(i % 20)
         if not res.is_valid:
-            color = "red"  # Highlight failing nets
+            color = "red"
 
         label_added = False
-        for _j, comp in enumerate(res.path):
-            # 1. Plot geometry
+        for comp in res.path:
+            # 1. Plot Collision Geometry (Translucent fill)
+            # This is the geometry used during search (e.g. proxy or arc)
             for poly in comp.geometry:
-                # Handle both Polygon and MultiPolygon (e.g. from SBend)
                 if isinstance(poly, MultiPolygon):
                     geoms = list(poly.geoms)
                 else:
                     geoms = [poly]
 
                 for g in geoms:
-                    x, y = g.exterior.xy
-                    ax.fill(x, y, alpha=0.7, fc=color, ec="black", label=net_id if not label_added else "")
+                    if hasattr(g, "exterior"):
+                        x, y = g.exterior.xy
+                        ax.fill(x, y, alpha=0.15, fc=color, ec=color, linestyle='--', lw=0.5, zorder=2)
+                    else:
+                        # Fallback for LineString or other geometries
+                        x, y = g.xy
+                        ax.plot(x, y, color=color, alpha=0.15, linestyle='--', lw=0.5, zorder=2)
+
+            # 2. Plot "Actual" Geometry (The high-fidelity shape used for fabrication)
+            # Use comp.actual_geometry if it exists (should be the arc)
+            actual_geoms_to_plot = comp.actual_geometry if comp.actual_geometry is not None else comp.geometry
+            
+            for poly in actual_geoms_to_plot:
+                if isinstance(poly, MultiPolygon):
+                    geoms = list(poly.geoms)
+                else:
+                    geoms = [poly]
+                for g in geoms:
+                    if hasattr(g, "exterior"):
+                        x, y = g.exterior.xy
+                        ax.plot(x, y, color=color, lw=1.5, alpha=0.9, zorder=3, label=net_id if not label_added else "")
+                    else:
+                        x, y = g.xy
+                        ax.plot(x, y, color=color, lw=1.5, alpha=0.9, zorder=3, label=net_id if not label_added else "")
                     label_added = True
 
-            # 2. Plot subtle port orientation arrow for internal ports
+            # 3. Plot subtle port orientation arrow
             p = comp.end_port
             rad = numpy.radians(p.orientation)
-            u = numpy.cos(rad)
-            v = numpy.sin(rad)
-            ax.quiver(p.x, p.y, u, v, color="black", scale=40, width=0.003, alpha=0.3, pivot="tail", zorder=4)
+            ax.quiver(p.x, p.y, numpy.cos(rad), numpy.sin(rad), color="black", 
+                      scale=40, width=0.002, alpha=0.2, pivot="tail", zorder=4)
 
-        # 3. Plot main arrows for netlist ports (if provided)
-        if netlist and net_id in netlist:
-            start_p, target_p = netlist[net_id]
+        if not res.path and not res.is_valid:
+            # Best-effort display: If the path is empty but failed, it might be unroutable.
+            # We don't have a partial path in RoutingResult currently.
+            pass
+
+    # 4. Plot main arrows for netlist ports
+    if netlist:
+        for net_id, (start_p, target_p) in netlist.items():
             for p in [start_p, target_p]:
                 rad = numpy.radians(p.orientation)
-                u = numpy.cos(rad)
-                v = numpy.sin(rad)
-                ax.quiver(p.x, p.y, u, v, color="black", scale=25, width=0.005, pivot="tail", zorder=6)
+                ax.quiver(p.x, p.y, numpy.cos(rad), numpy.sin(rad), color="black", 
+                          scale=25, width=0.004, pivot="tail", zorder=6)
 
     ax.set_xlim(bounds[0], bounds[2])
     ax.set_ylim(bounds[1], bounds[3])
     ax.set_aspect("equal")
-    ax.set_title("Inire Routing Results")
-    handles, labels = ax.get_legend_handles_labels()
-    if labels:
-        ax.legend()
-    plt.grid(True)
+    ax.set_title("Inire Routing Results (Dashed: Collision Proxy, Solid: Actual Geometry)")
+    
+    # Legend handling for many nets
+    if len(results) < 25:
+        handles, labels = ax.get_legend_handles_labels()
+        if labels:
+            ax.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize='small', ncol=2)
+            fig.tight_layout()
+
+    plt.grid(True, which='both', linestyle=':', alpha=0.5)
+    return fig, ax
+            
+def plot_danger_map(
+    danger_map: DangerMap,
+    ax: Axes | None = None,
+) -> tuple[Figure, Axes]:
+    """
+    Plot the pre-computed danger map as a heatmap.
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(10, 10))
+    else:
+        fig = ax.get_figure()
+
+    # Need to transpose because grid is [x, y] and imshow expects [row, col] (y, x)
+    # Also origin='lower' to match coordinates
+    im = ax.imshow(
+        danger_map.grid.T,
+        origin='lower',
+        extent=[danger_map.minx, danger_map.maxx, danger_map.miny, danger_map.maxy],
+        cmap='YlOrRd',
+        alpha=0.6
+    )
+    plt.colorbar(im, ax=ax, label='Danger Cost')
+    ax.set_title("Danger Map (Proximity Costs)")
+    return fig, ax
+
+def plot_expanded_nodes(
+    nodes: list[tuple[float, float, float]],
+    ax: Axes | None = None,
+    color: str = 'gray',
+    alpha: float = 0.3,
+) -> tuple[Figure, Axes]:
+    """
+    Plot A* expanded nodes for debugging.
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(10, 10))
+    else:
+        fig = ax.get_figure()
+
+    if not nodes:
+        return fig, ax
+
+    x, y, _ = zip(*nodes)
+    ax.scatter(x, y, s=1, c=color, alpha=alpha, zorder=0)
+    return fig, ax
+
+def plot_expansion_density(
+    nodes: list[tuple[float, float, float]],
+    bounds: tuple[float, float, float, float],
+    ax: Axes | None = None,
+    bins: int | tuple[int, int] = 50,
+) -> tuple[Figure, Axes]:
+    """
+    Plot a density heatmap (2D histogram) of expanded nodes.
+
+    Args:
+        nodes: List of (x, y, orientation) tuples.
+        bounds: (minx, miny, maxx, maxy) for the plot range.
+        ax: Optional existing axes to plot on.
+        bins: Number of bins for the histogram (int or (nx, ny)).
+
+    Returns:
+        Figure and Axes objects.
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(12, 12))
+    else:
+        fig = ax.get_figure()
+
+    if not nodes:
+        ax.text(0.5, 0.5, "No Expansion Data", ha='center', va='center', transform=ax.transAxes)
+        return fig, ax
+
+    x, y, _ = zip(*nodes)
+    
+    # Create 2D histogram
+    h, xedges, yedges = numpy.histogram2d(
+        x, y, 
+        bins=bins, 
+        range=[[bounds[0], bounds[2]], [bounds[1], bounds[3]]]
+    )
+    
+    # Plot as image
+    im = ax.imshow(
+        h.T, 
+        origin='lower', 
+        extent=[bounds[0], bounds[2], bounds[1], bounds[3]],
+        cmap='plasma',
+        interpolation='nearest',
+        alpha=0.7
+    )
+    
+    plt.colorbar(im, ax=ax, label='Expansion Count')
+    ax.set_title("Search Expansion Density")
+    ax.set_xlim(bounds[0], bounds[2])
+    ax.set_ylim(bounds[1], bounds[3])
+    
     return fig, ax