density.py

"""Windowed density map.

CMP planarization cares about LOCAL density, not just the layer total: a layer can hit its overall fill
target yet still have a sparse corner. This tiles the cell bounding box into a grid of windows and measures
the fill fraction of one layer in each window, by clipping the layer to the window and taking the covered
area over the window area. The result is a (ny, nx) array that densitydrc.py and report.py consume.
"""
from __future__ import annotations
import argparse
from dataclasses import dataclass

import numpy as np

import geom as _geom


@dataclass
class DensityMap:
    layer: tuple
    grid: np.ndarray          # (ny, nx) fill fraction per window, row 0 is the bottom row
    x_edges: np.ndarray       # length nx + 1
    y_edges: np.ndarray       # length ny + 1
    window_area: float        # area of one window (windows are equal-sized)


def _edges(lo, hi, n):
    return np.linspace(lo, hi, n + 1)


def _span(edges, lo, hi, n):
    """Inclusive [first, last] window index whose column/row overlaps the interval [lo, hi]."""
    i_lo = int(np.clip(np.searchsorted(edges, lo, side="right") - 1, 0, n - 1))
    i_hi = int(np.clip(np.searchsorted(edges, hi, side="right") - 1, 0, n - 1))
    return i_lo, i_hi


def _windows(bbox, nx, ny, window):
    xmin, ymin, xmax, ymax = bbox
    if window is not None:
        wx, wy = (window, window) if np.isscalar(window) else window
        nx = max(1, int(np.ceil((xmax - xmin) / wx)))
        ny = max(1, int(np.ceil((ymax - ymin) / wy)))
        return _edges(xmin, xmin + nx * wx, nx), _edges(ymin, ymin + ny * wy, ny)
    return _edges(xmin, xmax, nx), _edges(ymin, ymax, ny)


def density_map(layout, layer, nx=8, ny=8, window=None):
    """Fill fraction of `layer` in each window. Pass nx/ny window counts, or window=(wx, wy) in microns."""
    import gdstk

    layer = _geom._as_layer(layer)
    raw = layout.polygons.get(layer, [])
    x_edges, y_edges = _windows(layout.bbox, nx, ny, window)
    nx, ny = len(x_edges) - 1, len(y_edges) - 1
    grid = np.zeros((ny, nx), float)
    warea = (x_edges[1] - x_edges[0]) * (y_edges[1] - y_edges[0])
    if not raw or warea <= 0:
        return DensityMap(layer, grid, x_edges, y_edges, float(warea))

    # Bin each polygon into the windows its bounding box overlaps, so a window only clips the few
    # polygons that actually touch it rather than the whole layer.
    buckets = {}
    for pts in raw:
        a = np.asarray(pts, float)
        i_lo, i_hi = _span(x_edges, a[:, 0].min(), a[:, 0].max(), nx)
        j_lo, j_hi = _span(y_edges, a[:, 1].min(), a[:, 1].max(), ny)
        gp = gdstk.Polygon(a)
        for j in range(j_lo, j_hi + 1):
            for i in range(i_lo, i_hi + 1):
                buckets.setdefault((i, j), []).append(gp)

    for (i, j), polys in buckets.items():
        cell = gdstk.rectangle((x_edges[i], y_edges[j]), (x_edges[i + 1], y_edges[j + 1]))
        inter = gdstk.boolean(polys, [cell], "and")
        grid[j, i] = float(sum(p.area() for p in inter)) / warea
    return DensityMap(layer, grid, x_edges, y_edges, float(warea))


def _rect(x0, y0, x1, y1):
    return np.array([[x0, y0], [x1, y0], [x1, y1], [x0, y1]], float)


def _validate():
    # Left half of a 10x10 cell filled on layer (1,0): with 2x1 windows the left is full, the right empty.
    left = _geom.from_layers({(1, 0): [_rect(0, 0, 5, 10)]}, bbox=(0, 0, 10, 10))
    dm = density_map(left, (1, 0), nx=2, ny=1)
    ok_left = abs(dm.grid[0, 0] - 1.0) < 1e-9 and abs(dm.grid[0, 1] - 0.0) < 1e-9

    # Fully filled cell: every window reads 1.0.
    full = _geom.from_layers({(1, 0): [_rect(0, 0, 10, 10)]}, bbox=(0, 0, 10, 10))
    ok_full = np.allclose(density_map(full, (1, 0), nx=2, ny=2).grid, 1.0)

    # Bottom half filled, 1x2 windows: bottom window 1.0, top 0.0.
    bottom = _geom.from_layers({(1, 0): [_rect(0, 0, 10, 5)]}, bbox=(0, 0, 10, 10))
    db = density_map(bottom, (1, 0), nx=1, ny=2)
    ok_bottom = abs(db.grid[0, 0] - 1.0) < 1e-9 and abs(db.grid[1, 0] - 0.0) < 1e-9

    # A window that is half covered reads 0.5.
    half = _geom.from_layers({(1, 0): [_rect(0, 0, 10, 5)]}, bbox=(0, 0, 10, 10))
    ok_half = abs(density_map(half, (1, 0), nx=1, ny=1).grid[0, 0] - 0.5) < 1e-9

    for label, cond in [("left/right split", ok_left), ("uniform = 1.0", ok_full),
                        ("bottom/top split", ok_bottom), ("half window = 0.5", ok_half)]:
        print(f"[{label}]  {'PASS' if cond else 'FAIL'}")
    print("RESULT:", "PASS" if all([ok_left, ok_full, ok_bottom, ok_half]) else "FAIL")


def main():
    ap = argparse.ArgumentParser(description="windowed per-layer density map")
    ap.add_argument("--validate", action="store_true")
    ap.parse_args(); _validate()


if __name__ == "__main__":
    main()