import time
import warnings
from collections import defaultdict
from collections.abc import Sequence
from itertools import product
from typing import Literal, Protocol, TypedDict
import numpy as np
import tidy3d as td
from scipy.stats.qmc import LatinHypercube
from tidy3d.components.data.data_array import ScalarModeFieldDataArray
from .extension import (
Circle,
Component,
Interpolator,
MaskSpec,
Path,
Polygon,
Port,
PortSpec,
Rectangle,
Reference,
SMatrix,
Technology,
_content_repr,
_pack_rectangles,
boolean,
config,
frequency_classification,
grid_ceil,
grid_floor,
offset,
snap_to_grid,
)
from .typing import Frequency
_MonitorData = td.components.data.monitor_data.MonitorData
class _StatusDict(TypedDict, total=False):
progress: float
message: str
tasks: dict[str, object]
class _HasStatus(Protocol):
@property
def status(self) -> _StatusDict: ...
# Speed of light in vacuum (in µm/s)
C_0: float = 2.99792458e14
# Elementary charge (in C)
Q = 1.602176634e-19
# Planck's constant (in Js)
H = 6.62607015e-34
# Boltzmann constant (in J/K)
K_B = 1.380649e-23
# Number of points in an auto-derived fit grid, centered on the carrier. Odd so
# the grid samples the carrier itself (DC, for a baseband grid).
_DEFAULT_NUM_FREQUENCIES = 101
def _natural_port_key(name: str) -> tuple[str, int]:
prefix = name.rstrip("0123456789")
suffix = name[len(prefix) :]
return prefix, int(suffix) if suffix else -1
def _default_fit_frequencies(carrier_frequency: Frequency, time_step: float) -> np.ndarray:
"""Auto-derive a fit grid one Nyquist band wide, centered on the
carrier: ``[carrier - 0.5/dt, carrier + 0.5/dt]``. A baseband grid
(carrier 0) is the two-sided ``[-0.5/dt, 0.5/dt]``, handled like any
other carrier.
"""
nyquist = 0.5 / time_step
return carrier_frequency + np.linspace(-nyquist, nyquist, _DEFAULT_NUM_FREQUENCIES)
def _is_multiple_of_90(angle: float) -> bool:
rotation_fraction = angle % 90
return rotation_fraction < 1e-12 or (90 - rotation_fraction < 1e-12)
def _angles_equal(a: float, b: float) -> bool:
r = (a - b) % 360
return r <= 1e-12 or 360 - r <= 1e-12
def _gather_status(*runners: SMatrix | _HasStatus) -> _StatusDict:
"""Create an overall status based on a collection of runners."""
num_tasks = 0
progress = 0
message = "success"
tasks = {}
for task in runners:
task_status = (
{"progress": 100, "message": "success"} if isinstance(task, SMatrix) else task.status
)
inner_tasks = task_status.get("tasks", {})
tasks.update(inner_tasks)
task_weight = max(1, len(inner_tasks))
num_tasks += task_weight
if message != "error":
if task_status["message"] == "error":
message = "error"
elif task_status["message"] == "running":
message = "running"
progress += task_weight * task_status["progress"]
elif task_status["message"] == "success":
progress += task_weight * 100
if message == "running":
progress /= num_tasks
else:
progress = 100
return {"progress": progress, "message": message, "tasks": tasks}
def _align_and_overlap(
data0: _MonitorData, data1: _MonitorData, magnitude_warning: bool = True
) -> np.ndarray:
rotations = [(0, "+"), (1, "+"), (0, "-"), (1, "-")]
dir0 = getattr(data0.monitor, "direction", None)
if dir0 is None:
dir0 = data0.monitor.store_fields_direction
dir1 = getattr(data1.monitor, "direction", None)
if dir1 is None:
dir1 = data1.monitor.store_fields_direction
r0 = rotations.index((data0.monitor.size.index(0), dir0))
r1 = rotations.index((data1.monitor.size.index(0), dir1))
rotation = (r1 - r0) % 4
fields0 = data0.field_components
fields1 = data1.field_components
dims = fields0["Ez"].dims
coords = {d: fields0["Ez"].coords[d].values.copy() for d in dims}
center = (data0.monitor.center[0], data0.monitor.center[1])
if rotation == 0:
fields0 = {
"Ex": fields0["Ex"].values,
"Hx": fields0["Hx"].values,
"Ey": fields0["Ey"].values,
"Hy": fields0["Hy"].values,
"Ez": fields0["Ez"].values,
"Hz": fields0["Hz"].values,
}
elif rotation == 1:
x = coords["x"]
coords["x"] = -coords["y"]
coords["y"] = x
center = (-center[1], center[0])
ix = dims.index("x")
iy = dims.index("y")
fields0 = {
"Ex": np.swapaxes(-fields0["Ey"].values, ix, iy),
"Hx": np.swapaxes(-fields0["Hy"].values, ix, iy),
"Ey": np.swapaxes(fields0["Ex"].values, ix, iy),
"Hy": np.swapaxes(fields0["Hx"].values, ix, iy),
"Ez": np.swapaxes(fields0["Ez"].values, ix, iy),
"Hz": np.swapaxes(fields0["Hz"].values, ix, iy),
}
elif rotation == 2:
coords["x"] = -coords["x"]
coords["y"] = -coords["y"]
center = (-center[0], -center[1])
fields0 = {
"Ex": -fields0["Ex"].values,
"Hx": -fields0["Hx"].values,
"Ey": -fields0["Ey"].values,
"Hy": -fields0["Hy"].values,
"Ez": fields0["Ez"].values,
"Hz": fields0["Hz"].values,
}
elif rotation == 3:
x = coords["x"]
coords["x"] = coords["y"]
coords["y"] = -x
center = (center[1], -center[0])
ix = dims.index("x")
iy = dims.index("y")
fields0 = {
"Ex": np.swapaxes(fields0["Ey"].values, ix, iy),
"Hx": np.swapaxes(fields0["Hy"].values, ix, iy),
"Ey": np.swapaxes(-fields0["Ex"].values, ix, iy),
"Hy": np.swapaxes(-fields0["Hx"].values, ix, iy),
"Ez": np.swapaxes(fields0["Ez"].values, ix, iy),
"Hz": np.swapaxes(fields0["Hz"].values, ix, iy),
}
coords["x"] = coords["x"] + data1.monitor.center[0] - center[0]
coords["y"] = coords["y"] + data1.monitor.center[1] - center[1]
n, t = ("x", "y") if r1 % 2 == 0 else ("y", "x")
tangential_components = ("E" + t, "H" + t, "Ez", "Hz")
fields0 = {
c: ScalarModeFieldDataArray(fields0[c], dims=dims, coords=coords)
for c in tangential_components
}
coords1 = td.Coords(
x=fields1["Ez"].coords["x"].values,
y=fields1["Ez"].coords["y"].values,
z=fields1["Ez"].coords["z"].values,
)
fields0 = {c: coords1.spatial_interp(fields0[c], "linear") for c in tangential_components}
sign = -1 if r1 == 1 or r1 == 2 else 1
d_area = sign * data1._diff_area
e0_h1 = fields0["E" + t] * fields1["Hz"] - fields0["Ez"] * fields1["H" + t]
e1_h0 = fields1["E" + t] * fields0["Hz"] - fields1["Ez"] * fields0["H" + t]
integrand = (e0_h1 + e1_h0) * d_area
overlap = 0.25 * integrand.sum(dim=d_area.dims).isel({n: 0}, drop=True).values
# Modes are normalized by the mode solver, so the overlap should be only a phase difference.
# We normalize the result to remove numerical errors introduced by the grid interpolation.
overlap_mag = np.abs(overlap)
if magnitude_warning and not np.allclose(overlap_mag, 1.0, atol=0.1):
max_err = overlap_mag.flat[np.argmax(np.abs(overlap_mag - 1.0))]
warnings.warn(
f"Modal overlap calculation resulted in an unexpected magnitude ({max_err}). Consider "
"increasing the mesh refinement for the mode solver.",
RuntimeWarning,
2,
)
return overlap / overlap_mag
[docs]
def route_length(
component: Component, layer: Sequence[int] | None = None, port_spec: PortSpec | None = None
) -> float:
"""Measure the length of parametric routes.
Internally, this functions adds the path lengths (without offsets) for
all paths with distinct endpoints in a specific layer.
Args:
component: Component with routes to be measured.
layer: Layer used to look for paths. If ``None``, a best guess based
on ``port_spec`` will be used.
port_spec: Port specification used for the route. If ``None``, the
component will be inspected and a best guess used.
Returns:
Total path length.
See also:
- `Parametric routes <../parametric.rst#routing>`__
- :func:`effective_route_length`
"""
if layer is None:
if port_spec is None:
ports = tuple(p for p in component.ports.values() if isinstance(p, Port))
opt_ports = tuple(p for p in ports if p.classification == "optical")
elec_ports = tuple(p for p in ports if p.classification == "electrical")
if len(opt_ports) == 2:
port_spec = opt_ports[0].spec
elif len(elec_ports) == 2:
port_spec = elec_ports[0].spec
if port_spec is not None:
profiles = sorted(port_spec.path_profiles_list(), key=lambda p: (abs(p[1]), p[0], p[2]))
if len(profiles) > 0:
layer = profiles[0][2]
structures = component.get_structures(layer)
if layer is not None:
structures = {(0, 0): structures}
result = 0.0
for structure_list in structures.values():
paths = defaultdict(list)
for path in structure_list:
if isinstance(path, Path):
key = tuple(
sorted(
tuple(snap_to_grid(pos))
for pos in (path.origin, path.at(path.size, output="position"))
)
)
paths[key].append(path)
if len(paths) > 0:
result = max(
result,
sum(
sum(path.length(include_offset=False) for path in path_list) / len(path_list)
for path_list in paths.values()
),
)
return result
def _log_samples(a: float, b: float, r: float) -> np.ndarray:
q = b / a
n = 1 + np.ceil(np.log(q) / np.log(1 + r))
if n < 2:
return np.array([a, b]) if a != b else np.array([a])
r = q ** (1.0 / (n - 1))
return a * r ** np.arange(n)
def _samples_2d(interp: Interpolator) -> np.ndarray:
x = interp.x
return x if x.ndim == 2 else x.reshape(-1, 1)
[docs]
def effective_route_length(
*,
paths: Sequence[Path] = (),
component: Component | None = None,
layer: Sequence[int] | None = None,
port_spec: PortSpec | None = None,
n_eff: Interpolator | None = None,
frequency: Frequency | None = None,
technology: Technology | None = None,
width_rtol: float = 0.2,
offset_rtol: float = 0.2,
curvature_rtol: float = 0.2,
verbose: bool = True,
show_progress: bool = True,
) -> float:
r"""Measure the effective length of waveguides.
The effective length of a waveguide is the integral of its effective
index along the propagation length:
.. math:: \int n_\text{eff}(s) \, {\rm d}s
If provided, the ``n_eff`` interpolator can be parametrized by
"frequency" or "wavelength", bend "radius" or "curvature", waveguide
"width", "offset", and "angle" (in degrees, useful in anisotropic
media). The values used for sampling and numerical integration are
derived from the ``paths`` provided to the function or found in the
``component``.
When ``n_eff`` is not provided, the paths are inspected to define the
parameter ranges for "width", "curvature", and "angle" (if anisotropic
media is used in ``technology``). Then mode solver runs based on
``port_spec`` are used to evaluate the effective index and create an
interpolator. In particular, the "width" variation is only applied to
a single path profile in ``port_spec``, if one can be detected.
Args:
paths: Paths to use for effective length calculation.
component: Component to look for paths if ``paths`` is not given.
layer: Layer used to look for paths. If ``None``, a best guess based
on ``port_spec`` will be used.
port_spec: Port specification used for the base waveguide. If
``None``, the component will be inspected and a best guess used.
n_eff: Interpolator with parametrized effective index values.
frequency: Frequency for effective index computation. It must be
provided if ``n_eff`` is not provided, or if it is parametrized by
"frequency" or "wavelength".
technology: Technology used for mode-solving, if needed.
width_rtol: Maximum allowed relative change in path width between
consecutive integration points. If less than or equal to 0, width
parametrization is disabled.
offset_rtol: Maximum allowed relative change in offset between
consecutive integration points. If less than or equal to 0, offset
parametrization is disabled.
curvature_rtol: Maximum allowed relative change in curvature between
consecutive integration points. If less than or equal to 0,
curvature parametrization is disabled.
verbose: Flag controlling mode solver verbosity.
show_progress: Flag controlling mode solver progress.
Returns:
Effective path length for guided waves.
See also:
:func:`route_length`
"""
requested_layer = layer
if technology is None:
technology = config.default_technology if component is None else component.technology
if isinstance(port_spec, str):
port_spec = technology.ports[port_spec]
if component is not None and (len(paths) == 0 or n_eff is None):
if port_spec is None:
ports = (
component.ports
if frequency is None
else component.select_ports(frequency_classification(frequency))
)
opt_ports = tuple(
p for p in ports.values() if isinstance(p, Port) and p.classification == "optical"
)
elec_ports = tuple(
p
for p in ports.values()
if isinstance(p, Port) and p.classification == "electrical"
)
if layer is not None:
opt_ports = tuple(
p
for p in opt_ports
if any(
(pp[0] < p.spec.width and pp[1] == 0 and pp[2] == layer)
for pp in p.spec.path_profiles_list()
)
)
elec_ports = tuple(
p
for p in elec_ports
if any(
(pp[0] < p.spec.width and pp[1] == 0 and pp[2] == layer)
for pp in p.spec.path_profiles_list()
)
)
if len(opt_ports) == 2 and opt_ports[0].can_connect_to(opt_ports[1]):
port_spec = opt_ports[0].spec
elif len(elec_ports) == 2 and elec_ports[0].can_connect_to(elec_ports[1]):
port_spec = elec_ports[0].spec
if port_spec is not None and layer is None:
profiles = sorted(
(p for p in port_spec.path_profiles_list() if p[0] < port_spec.width),
key=lambda p: (abs(p[1]), -p[0], p[2]),
)
if len(profiles) > 0:
layer = profiles[0][2]
if layer is None:
all_structures = component.get_structures()
if len(all_structures) > 0:
layer = max(
all_structures,
key=lambda k: sum(s.length() for s in all_structures[k] if isinstance(s, Path)),
)
if len(paths) == 0 and component is not None:
if layer is None:
# This can only happen if len(component.get_structures()) == 0
return 0.0
paths = {}
for path in component.get_structures(layer):
if isinstance(path, Path):
key = tuple(
sorted(
tuple(snap_to_grid(pos))
for pos in (path.origin, path.at(path.size, output="position"))
)
)
paths[key] = path
paths = list(paths.values())
if n_eff is None:
from .models.tidy3d import ( # noqa: PLC0415
_isotropic_uniform,
_local_tidy3d_solver,
_ModeSolverRunner,
)
if port_spec is None or frequency is None:
raise RuntimeError("Please provide either 'n_eff', or 'port_spec' and 'frequency'.")
classification = frequency_classification(frequency)
clad_profiles = []
core_profiles = []
for p in port_spec.path_profiles_list():
if (
p[1] == 0.0
and p[0] < port_spec.width
and (requested_layer is None or p[2] == requested_layer)
):
core_profiles.append(list(p))
else:
clad_profiles.append(list(p))
adjust_width = len(core_profiles) == 1
path_profiles = core_profiles + clad_profiles
width_min = np.inf
width_max = -np.inf
curv_min = np.inf
curv_max = 0
for path in paths:
s, _, _, curvature, width, _ = path._integration_points(
width_rtol=width_rtol,
offset_rtol=offset_rtol,
curvature_rtol=curvature_rtol,
include_offset=False,
)
if s.size < 2:
continue
w_min = width.min()
if w_min < width_min:
width_min = w_min
w_max = width.max()
if w_max > width_max:
width_max = w_max
curvature = np.abs(curvature[curvature != 0.0])
if curvature.size > 0:
c_min = curvature.min()
if c_min < curv_min:
curv_min = c_min
c_max = curvature.max()
if c_max > curv_max:
curv_max = c_max
if not adjust_width and width_max > width_min:
raise RuntimeError(
"Unable to identify waveguide core among the path profiles of the port "
"specification to generate width variations. Please provide 'n_eff' directly to "
"include width variations."
)
names = []
ranges = []
if width_rtol > 0 and width_max > width_min:
names.append("width")
w_min = grid_floor(width_min)
# use grid if w_min == 0
ranges.append(
_log_samples(max(config.grid, w_min), grid_ceil(width_max), 3 * width_rtol)
)
if w_min < config.grid:
ranges[-1] = np.hstack(([w_min], ranges[-1]))
elif adjust_width:
path_profiles[0][0] = width_min
if curvature_rtol > 0 and curv_max > 0:
names.append("curvature")
curv_samples = _log_samples(curv_min * 0.99, curv_max * 1.01, 3 * curvature_rtol)
ranges.append(np.hstack(([0.0], curv_samples)))
isotropic = _isotropic_uniform(technology, classification)
if not isotropic:
names.append("angle")
ranges.append(np.linspace(0, 90, 7))
ndim = len(ranges)
if ndim == 0:
names = ["curvature"]
samples = [[0.0]]
elif ndim == 1:
samples = ranges[0].reshape((-1, 1))
else:
lengths = [len(s) for s in ranges]
corners = 2**ndim
random = max(1, max(lengths) - ndim)
samples = np.empty((corners + random, ndim), dtype=float)
samples[:corners, :] = list(product(*((s[0], s[-1]) for s in ranges)))
x = LatinHypercube(ndim, rng=np.random.default_rng(1)).random(random)
for i, (name, r) in enumerate(zip(names, ranges, strict=True)):
if name == "width":
samples[corners:, i] = r[0] ** (1 - x[:, i]) * r[-1] ** x[:, i]
elif name == "curvature":
samples[corners:, i] = r[1] ** (1 - x[:, i]) * r[-1] ** x[:, i]
else:
samples[corners:, i] = r[0] * (1 - x[:, i]) + r[-1] * x[:, i]
if show_progress:
print("Starting…", end="\r", flush=True)
runners = []
port_spec = port_spec.copy()
port_spec.path_profiles = path_profiles
port = Port((0, 0), 0, port_spec)
for sample in samples:
for i, name in enumerate(names):
if name == "width":
path_profiles[0][0] = sample[i]
port_spec.path_profiles = path_profiles
elif name == "curvature":
port.bend_radius = 0 if sample[i] == 0.0 else 1 / sample[i]
elif name == "angle":
port.input_direction = sample[i]
mode_solver = port.to_tidy3d_mode_solver(
[frequency], None, False, technology, isotropic
)
runners.append(
_ModeSolverRunner(mode_solver, [frequency], None, technology, verbose=verbose)
)
# Wait for all to finish
progress_chars = "-\\|/"
i = 0
while True:
status = _gather_status(*runners)
message = status["message"]
if message == "success":
if show_progress:
print("Progress: 100% ", end="\n", flush=True)
break
elif message == "running":
if _local_tidy3d_solver:
raise RuntimeError("Unexpected on-prem mode solver message.")
if show_progress:
p = max(0, min(100, int(status.get("progress", 0))))
c = progress_chars[i]
i = (i + 1) % len(progress_chars)
print(f"Progress: {p}% {c}", end="\r", flush=True)
time.sleep(0.3)
elif message == "error":
if show_progress:
print("Progress: error", end="\n", flush=True)
raise RuntimeError("Mode solver runs resulted in error.")
else:
raise RuntimeError(f"Status message unknown: {message!r}.")
values = [r.data.n_eff.values[0, 0] for r in runners]
n_eff = Interpolator(samples, values, parameter_names=names)
if np.ndim(n_eff.y) != 1:
raise TypeError("'n_eff' interpolator must have a single value.")
if n_eff.parameter_names is None:
raise RuntimeError(
"'n_eff' interpolator must have parameter names to correctly map path properties."
)
parameters = {n: i for i, n in enumerate(n_eff.parameter_names)}
# Replace wavelength with frequency to match argument
if "wavelength" in parameters:
i = parameters["wavelength"]
samples = _samples_2d(n_eff)
if "frequency" in parameters:
samples = np.hstack((samples[:, :i], samples[:, i + 1 :]))
names = list(n_eff.parameter_names)
names.remove("wavelength")
n_eff = Interpolator(
samples, n_eff.y, n_eff.method, n_eff.coords, n_eff.extrapolation, names
)
else:
samples[:, i] = C_0 / samples[:, i]
names = list(n_eff.parameter_names)
names[i] = "frequency"
n_eff = Interpolator(
samples, n_eff.y, n_eff.method, n_eff.coords, n_eff.extrapolation, names
)
parameters = {n: i for i, n in enumerate(n_eff.parameter_names)}
# Replace radius with curvature for correct interpolation of R → ∞
if "radius" in parameters:
i = parameters["radius"]
samples = _samples_2d(n_eff)
if "curvature" in parameters:
samples = np.hstack((samples[:, :i], samples[:, i + 1 :]))
names = list(n_eff.parameter_names)
names.remove("radius")
n_eff = Interpolator(
samples, n_eff.y, n_eff.method, n_eff.coords, n_eff.extrapolation, names
)
else:
r = samples[:, i]
r[r == 0.0] = np.inf
samples[:, i] = 1.0 / r
names = list(n_eff.parameter_names)
names[i] = "curvature"
n_eff = Interpolator(
samples, n_eff.y, n_eff.method, n_eff.coords, n_eff.extrapolation, names
)
parameters = {n: i for i, n in enumerate(n_eff.parameter_names)}
for name in parameters:
if name not in ("width", "offset", "angle", "curvature", "frequency"):
raise RuntimeError(
f"'n_eff' interpolator parameter {name!r} is not supported in this function."
)
width_index = parameters.get("width")
if width_index is None:
width_rtol = 0
offset_index = parameters.get("offset")
if offset_index is None:
offset_rtol = 0
curvature_index = parameters.get("curvature")
if curvature_index is None:
curvature_rtol = 0
else:
# If anisotropic, curvature sign matters, otherwise, it shouldn't. Try to guess user intent.
absolute_curvature = np.all(_samples_2d(n_eff)[:, curvature_index] >= 0.0)
angle_index = parameters.get("angle")
if angle_index is not None:
samples = _samples_2d(n_eff)
angles = samples[:, angle_index] % 360
samples[:, angle_index] = angles
angle_max = angles.max()
angle_min = angles.min()
angle_range = None
if angle_min >= 0:
if angle_max <= 90:
angle_range = 90
elif angle_max <= 180:
angle_range = 180
values = n_eff.y
if angle_range is None:
mask = angles == 0.0
if mask.any():
new_samples = np.array(samples[mask])
new_samples[:, angle_index] = 360.0
samples = np.vstack((samples, new_samples))
values = np.hstack((values, values[mask]))
n_eff = Interpolator(
samples, values, n_eff.method, n_eff.coords, n_eff.extrapolation, n_eff.parameter_names
)
frequency_index = parameters.get("frequency")
if frequency_index is not None and frequency is None:
raise RuntimeError("Interpolator requires 'frequency' argument.")
result = 0.0
for path in paths:
s, _, direction, curvature, width, offset = path._integration_points(
width_rtol=width_rtol,
offset_rtol=offset_rtol,
curvature_rtol=curvature_rtol,
include_offset=False,
)
if s.size < 2:
continue
x = np.empty((s.size, len(parameters)), dtype=float)
if frequency_index is not None:
x[:, frequency_index] = frequency
if width_index is not None:
x[:, width_index] = width
if offset_index is not None:
x[:, offset_index] = offset
if curvature_index is not None:
x[:, curvature_index] = np.abs(curvature) if absolute_curvature else curvature
if angle_index is not None:
angle = np.arctan2(direction[:, 1], direction[:, 0]) / np.pi * 180.0 # [-180; 180]
if angle_range is not None:
angle = np.abs(angle) # [0; 180]
if angle_range == 90:
angle[angle > 90] = 180 - angle[angle > 90] # [0; 90]
else:
angle = angle % 360 # [0; 360)
x[:, angle_index] = angle
y = n_eff(x if len(n_eff.parameter_names) > 1 else x[:, 0])
ds = np.diff(s)
result += (0.5 * (y[:-1] + y[1:]) * ds).sum()
return result
def _detect_routes(component):
match component.parametric_function:
case "photonforge.parametric.route" | "photonforge.parametric.route_s_bend":
return [component]
case (
"photonforge.parametric.route_u"
| "photonforge.parametric.route_z"
| "photonforge.parametric.route_l"
| "photonforge.parametric.route_auto"
):
return (
[component]
if len(component.ports) == 2
else [ref.component for ref in component.references]
)
case _:
warnings.warn(
f"Component {component.name!r} is not recognized as a route created from the "
f"internal parametric library. 'routing_collisions' may not work as expected.",
RuntimeWarning,
3,
)
return [component]
def _collect_obstacles(obj, collision_layers):
if isinstance(obj, (Rectangle, Circle, Path, Polygon)):
yield obj.to_polygon()
elif isinstance(obj, (Component, Reference)):
if collision_layers is not None:
for layer in collision_layers:
yield from (s.to_polygon() for s in obj.get_structures(layer))
else:
for structs in obj.get_structures().values():
yield from (s.to_polygon() for s in structs)
else:
for inner in obj:
yield from _collect_obstacles(inner, collision_layers)
def _intersects(a, b):
return not (a[1][0] < b[0][0] or b[1][0] < a[0][0] or a[1][1] < b[0][1] or b[1][1] < a[0][1])
[docs]
def routing_collisions(route, obstacles=(), collision_layers=None, collision_offset=0):
"""Find route collisions with obstacles and sibling route segments.
The route should be a component created by the internal parametric
routing functions. Obstacles may be 2D structures, components,
references, or sequences of those objects. Components and references are
expanded through :meth:`get_structures`, optionally limited to
``collision_layers``.
Args:
route: Route component to check.
obstacles: Obstacles to test against the route geometry.
collision_layers: Layers used to collect route and obstacle
geometry. If ``None``, all structure layers are used.
collision_offset: Offset applied to route geometry before testing.
Returns:
Tuple ``(obstacle_intersections, route_crossings)``. The first list
contains route-obstacle intersection polygons. The second list
contains intersections between route segments.
"""
routes = [
[(p.bounds(), p) for p in _collect_obstacles(r, collision_layers)]
for r in _detect_routes(route)
]
all_obstacles = [(p.bounds(), p) for p in _collect_obstacles(obstacles, collision_layers)]
obstacle_intersections = []
route_crossings = []
for i, route_data in enumerate(routes):
offset_route = offset([p for _, p in route_data], collision_offset)
offset_bounds = [p.bounds() for p in offset_route]
obstacles = [p for b, p in all_obstacles if any(_intersects(a, b) for a in offset_bounds)]
obstacle_intersections.extend(boolean(offset_route, obstacles, "*"))
for j in range(i):
route = [r for b, r in routes[j] if any(_intersects(a, b) for a in offset_bounds)]
if len(route) > 0:
route_crossings.extend(boolean(offset_route, route, "*"))
return obstacle_intersections, route_crossings
def _layer_in_mask_score(layer: tuple[int, int], mask: MaskSpec) -> int:
if mask.layer is not None:
return 1 if mask.layer == layer else None
operands = mask.operand1 + mask.operand2
if mask.operation == "+":
for inner in operands:
score = _layer_in_mask_score(layer, inner)
if score is not None:
return 10 * score + len(operands)
elif mask.operation == "*":
for inner in operands:
score = _layer_in_mask_score(layer, inner)
if score is not None:
return 20 * score + len(operands)
elif mask.operation == "-":
for inner in mask.operand1:
score = _layer_in_mask_score(layer, inner)
if score is not None:
return 20 * score + len(operands)
return None
_virtual_port_specs = {}
[docs]
def virtual_port_spec(
num_modes: int = 1, classification: str = "optical", impedance: complex | Interpolator = 50
) -> PortSpec:
"""Template to generate a virtual PortSpec.
Virtual port specs have no path profiles and can be used to help with
schematic-driven design before any layout is created.
Args:
num_modes: Number of modes supported by the port.
classification: One of ``"optical"`` or ``"electrical"``.
impedance: Complex impedance as a single or frequency-dependent
interpolated value (in ohms).
Returns:
Virtual port specification with no path profiles.
"""
virtual = None
if classification == "optical":
virtual = PortSpec("Virtual spec (optical)", 1, (0, 0), num_modes)
elif classification == "electrical":
virtual = PortSpec("Virtual spec (electrical)", 1, (0, 0), num_modes, impedance=impedance)
key = _content_repr(classification, num_modes, impedance, include_config=False)
cached = _virtual_port_specs.get(key)
if cached != virtual:
_virtual_port_specs[key] = virtual
cached = virtual
return cached
[docs]
def cpw_spec(
layer: str | Sequence[int],
signal_width: float,
gap: float,
ground_width: float | None = None,
description: str | None = None,
width: float | None = None,
limits: Sequence[float] | None = None,
num_modes: int = 1,
added_solver_modes: int = 0,
target_neff: float = 4.0,
gap_layer: None | str | Sequence[int] | None = None,
include_ground: bool = True,
conductor_limits: Sequence[float] | None = None,
technology: Technology | None = None,
) -> PortSpec:
"""Template to generate a coplanar transmission line PortSpec.
Args:
layer: Layer used for the transmission line layout.
signal_width: Width of the central conductor.
gap: Distance between the central conductor and the grounds.
ground_width: Width of the ground conductors.
description: Description used in :attr:`PortSpec.description`.
width: Dimension used in :attr:`PortSpec.width`.
limits: Vertical port limits used in :attr:`PortSpec.limits`.
num_modes: Value used for :attr:`PortSpec.num_modes`.
added_solver_modes: Value used for
:attr:`PortSpec.added_solver_modes`.
target_neff: Value used for :attr:`PortSpec.target_neff`.
gap_layer: If set, path profiles for the gap region are included in
this layer.
include_ground: If ``False``, ground path profiles are not included.
conductor_limits: Lower and upper bounds of the conductor layer
extrusion.
technology: Technology in use. If ``None``, the default is used.
Returns:
PortSpec for the CPW transmission line.
Note:
If ``conductor_limits`` is not given, the extrusion specifications
in ``technology`` are inspected. If an specification for the
selected ``layer`` is found, its extrusion limits are used.
"""
if technology is None:
technology = config.default_technology
if isinstance(layer, str):
layer = technology.layers[layer].layer
if isinstance(gap_layer, str):
gap_layer = technology.layers[gap_layer].layer
if conductor_limits is None:
best_score = 1e30
for extrusion in technology.extrusion_specs:
medium = extrusion.get_medium("electrical")
if not (medium.is_pec or isinstance(medium, td.LossyMetalMedium)):
continue
if extrusion.mask_spec.layer == layer:
conductor_limits = extrusion.limits
break
score = _layer_in_mask_score(layer, extrusion.mask_spec)
if score is not None and score < best_score:
conductor_limits = extrusion.limits
best_score = score
if conductor_limits is None:
raise RuntimeError(
f"Unable to find a conductor extrusion specification for layer {layer}. Please "
f"specify 'conductor_limits' manually."
)
z_center = 0.5 * (conductor_limits[0] + conductor_limits[1])
z_thickness = abs(conductor_limits[1] - conductor_limits[0])
cpw_min = min(signal_width, gap, z_thickness)
# Scale found manually by testing a range of configurations
cpw_scale = gap**0.3 * signal_width**0.6
ground_factor = 10
z_factor = 12
if ground_width is None:
ground_width = ground_factor * cpw_scale
offset = (signal_width + ground_width) / 2 + gap
full_width = signal_width + 2 * gap + 2 * ground_width
if description is None:
description = f"CPW (signal width: {signal_width}, gap: {gap})"
if width is None:
width = min(full_width, signal_width + 2 * (gap + ground_factor * cpw_scale)) - cpw_min
elif width >= full_width:
warnings.warn(
"CPW width is larger than the ground conductor extension. Please increase "
"'ground_width' or decrease 'width', otherwise check the port modes to "
"make sure the mode solver finds the correct modes.",
stacklevel=2,
)
if limits is None:
z_margin = z_thickness / 2 + z_factor * cpw_scale
limits = (z_center - z_margin, z_center + z_margin)
path_profiles = {"signal": (signal_width, 0, layer)}
if include_ground:
path_profiles["gnd0"] = (ground_width, -offset, layer)
path_profiles["gnd1"] = (ground_width, offset, layer)
if gap_layer is not None:
gap_offset = (signal_width + gap) / 2
path_profiles["gap0"] = (gap, -gap_offset, gap_layer)
path_profiles["gap1"] = (gap, gap_offset, gap_layer)
return PortSpec(
description=description,
width=width,
limits=limits,
num_modes=num_modes,
added_solver_modes=added_solver_modes,
target_neff=target_neff,
path_profiles=path_profiles,
voltage_path=[(signal_width / 2 + gap, z_center), (signal_width / 2, z_center)],
current_path=Rectangle(center=(0, z_center), size=(signal_width + gap, z_thickness + gap)),
)
[docs]
def grid_layout(
objects: Sequence[Component | Reference | Circle | Path | Polygon | Rectangle],
gap: float | Sequence[float] = 0,
shape: Sequence[int] | None = None,
align_x: Literal["left", "right", "center", "origin"] | None = "center",
align_y: Literal["bottom", "top", "center", "origin"] | None = "center",
direction: Literal[
"lr-bt", "lr-tb", "rl-bt", "rl-tb", "bt-lr", "tb-lr", "bt-rl", "tb-rl"
] = "lr-bt",
include_ports: bool = True,
layer: tuple[int] = (0, 0),
name: str | None = None,
) -> Component:
"""
Arrange components or other structures in a grid layout.
Args:
objects: Sequence of objects to arrange. They can be instances of
:class:`Component`, :class:`Reference`, or 2D structures.
gap: Horizontal and vertical gaps added between objects.
shape: Grid shape, specified as ``(columns, rows)``.
align_x: Horizontal alignment within the grid cell.
align_y: Vertical alignment within the grid cell.
direction: Placement order in the grid. Must be a combination of
``"lr"`` (left-to-right) or ``"rl"`` (right-to-left), and
``"bt"`` (bottom-to-top) or ``"tb"`` (top-to-bottom), as in
``"lr-bt"``, ``"rl-tb"``, ``"tb-lr"``, etc.
include_ports: Whether or not to include ports when computing
component bounds.
layer: If arraging geometrical structures, add them to this layer.
name: Name of the resulting component.
Returns:
Component with the objects arranged in a grid.
"""
num_objects = len(objects)
if num_objects == 0:
raise RuntimeError("List of objects cannot be empty.")
directions = {"rl-bt", "rl-tb", "lr-bt", "lr-tb", "bt-rl", "tb-rl", "bt-lr", "tb-lr"}
if direction not in directions:
alternatives = ", ".join(repr(d) for d in sorted(directions))
raise ValueError(f"Invalid value for 'direction'. Must be one of {alternatives}")
rows = int(num_objects**0.5 + 0.5) if shape is None else shape[1]
cols = (num_objects + rows - 1) // rows if shape is None else shape[0]
if num_objects > rows * cols:
raise ValueError("More components than available grid slots.")
if name is None:
name = f"GRID_{cols}_{rows}"
bounds = np.array(
[
obj.bounds(include_ports) if isinstance(obj, Component) else obj.bounds()
for obj in objects
]
)
size = (bounds[:, 1, :] - bounds[:, 0, :]).max(axis=0)
if align_x == "origin":
size[0] = bounds[:, 1, 0].max() - bounds[:, 0, 0].min()
if align_y == "origin":
size[1] = bounds[:, 1, 1].max() - bounds[:, 0, 1].min()
size += gap
x_offsets = [col * size[0] for col in range(cols)]
y_offsets = [row * size[1] for row in range(rows)]
if direction[:2] == "rl" or direction[3:] == "rl":
x_offsets.reverse()
if direction[:2] == "tb" or direction[3:] == "tb":
y_offsets.reverse()
if "r" in direction[:2]:
offsets = ((x, y) for y in y_offsets for x in x_offsets)
else:
offsets = ((x, y) for x in x_offsets for y in y_offsets)
offsets = np.array(tuple(offsets)[:num_objects])
if align_x == "left":
offsets[:, 0] -= bounds[:, 0, 0]
elif align_x == "right":
offsets[:, 0] -= bounds[:, 1, 0]
elif align_x == "center":
offsets[:, 0] -= bounds[:, :, 0].sum(axis=1) / 2
if align_y == "bottom":
offsets[:, 1] -= bounds[:, 0, 1]
elif align_y == "top":
offsets[:, 1] -= bounds[:, 1, 1]
elif align_y == "center":
offsets[:, 1] -= bounds[:, :, 1].sum(axis=1) / 2
technology = None
for i in range(len(objects)):
if isinstance(objects[i], Component):
objects[i] = Reference(objects[i])
if technology is None and isinstance(objects[i], Reference):
technology = objects[i].component.technology
objects[i].translate(offsets[i])
c = Component(name, technology)
c.add(layer, *objects)
return c
[docs]
def pack_layout(
objects: Sequence[Component | Reference | Circle | Path | Polygon | Rectangle],
gap: float | Sequence[float] = 0,
max_size: Sequence[float] = (0, 0),
aspect_ratio: float = 0,
grow_factor: float = 1.1,
sorting: Literal["best", "area"] | None = "best",
allow_rotation: bool = False,
method: Literal["bl", "blsf", "bssf", "baf", "cp"] = "blsf",
include_ports: bool = True,
layer: tuple[int] = (0, 0),
name: str = "PACK_{i}",
) -> list[Component]:
"""
Arrange components or other structures in a grid layout.
Args:
objects: Sequence of objects to arrange. They can be instances of
:class:`Component`, :class:`Reference`, or 2D structures.
gap: Horizontal and vertical gaps added between objects.
max_size: Maximal size of the packed component. If not all objects
fit in a single pack, multiple are used.
aspect_ratio: Desired width:height ratio for the pack.
grow_factor: Controls pack size increment. Values closer to 1 can
result in tighter packs at the cost of more computation.
sorting: Sorting option for the list of objects. If ``None``,
objects are packed in the order they are listed; ``'area'`` will
pack from largest to smallest, and ``"best"`` will try to choose
the best object to pack at each iteration.
allow_rotation: If ``True``, objects may be rotated by 90°.
method: Heuristic used to select a free slot during packing. See
below for information about the options.
include_ports: Whether or not to include ports when computing
component bounds.
layer: If arranging geometrical structures, add them to this layer.
name: Name template for the resulting components. Variable ``i`` is
used to indicate the pack index in the case of multiple packs.
Returns:
List of components with the packed objects.
Note:
The available methods for selecting a free slot for packing are:
Bottom left rule (``"bl"``):
Use the left-most position among the lowest upper y-value options.
Best long side fit (``"blsf"``):
Use the position that minimizes the leftover length on the long
side.
Best short side fit (``"bssf"``):
Use the position that minimizes the leftover length on the short
side.
Best area fit (``"baf"``):
Use the smallest available area that fits.
Contact point rule (``"cp"``):
Use the position that maximizes the length of the perimeter that
touches other objects.
Reference: Jukka Jylänki, *A Thousand Ways to Pack the Bin – A Practical
Approach to Two-Dimensional Rectangle Bin Packing*, 2010.
"""
if len(objects) == 0:
raise RuntimeError("List of objects cannot be empty.")
technology = None
for i in range(len(objects)):
if isinstance(objects[i], Component):
technology = objects[i].technology
break
if isinstance(objects[i], Reference):
technology = objects[i].component.technology
break
bounds = np.array(
[
obj.bounds(include_ports) if isinstance(obj, Component) else obj.bounds()
for obj in objects
]
)
sizes = bounds[:, 1, :] - bounds[:, 0, :] + gap
keep_order = sorting != "best"
if sorting == "area":
order = sorted(((a * b, i) for i, (a, b) in enumerate(sizes)), reverse=True)
objects = [objects[i] for _, i in order]
sizes = [sizes[i] for _, i in order]
else:
objects = list(objects)
sizes = list(sizes)
for i in range(2):
if max_size[i] > 0 and any(size[i] > max_size[i] for size in sizes):
for j in range(len(objects)):
if sizes[j][i] > max_size[i]:
raise RuntimeError(
f"{('Width', 'Height')[i]} of 'objects[{j}]' (plus gap) is larger than "
f"'max_size[{i}]' ({sizes[j][i]} > {max_size[i]})."
)
packs = []
while len(objects) > 0:
pack = _pack_rectangles(
sizes, method, max_size, aspect_ratio, grow_factor, keep_order, allow_rotation
)
if len(pack) == 0:
raise RuntimeError("Unable to pack objects.")
packed_objects = []
for index, corner, rotate in pack:
obj = objects[index]
objects[index] = None
sizes[index] = None
if isinstance(obj, Component):
obj = Reference(obj)
if rotate:
obj.rotate(90)
xy_min, _ = obj.bounds()
obj.translate(corner - xy_min)
packed_objects.append(obj)
c = Component(name.format(i=len(packs)), technology)
c.add(layer, *packed_objects)
packs.append(c)
objects = [obj for obj in objects if obj is not None]
sizes = [size for size in sizes if size is not None]
return packs