"""Boundary-conforming Delaunay layout inspired by DistMesh.
Reference: https://doi.org/10.1137/S0036144503429121
"""
from __future__ import annotations
import numpy as np
import shapely
from .._validation import _validate_integer, _validate_positive
from ..boundaries.base import _inset_interior, _project_inside
from .base import (
LayoutAlgorithm,
_allocate_counts,
_kept_apart_mask,
_nearest_neighbor_distances,
_spread_selection,
)
from .perimeter import _contained_perimeter_points
__all__ = ["Delaunay"]
#: Share of the arc-length interval below which a derived boundary node is
#: too close to an earlier one and becomes an interior node.
_CROWDED_BOUNDARY = 0.8
#: Share of the hexagonal packing count a fill's trial mesh draws. A settled
#: mesh keeps its closest pair near 0.94 of the target spacing, so 0.85 to
#: 0.89 of the packing count fits over offshore shapes, and the trial lands
#: near the answer.
_TRIAL_SHARE = 0.85
#: Share of the median nearest-neighbor distance a settled mesh keeps as its
#: closest pair at the count a fill holds. Over 80 minimum-spacing fills of
#: offshore and other shapes 0.996 predicts the count found within half a
#: turbine on average, and the search needs 2.6 meshes after the trial,
#: where 0.99 fell four counts short and needed 6.
_CLOSEST_SHARE = 0.996
#: Share of the target spacing below which an edge repels its nodes. The
#: target is the side of ``n_turbines`` equilateral triangles over the area.
#: A fixed length lets the mesh settle once no edge is shorter. A length
#: scaled with the RMS edge length, 1.2 times it in DistMesh and 0.85 times
#: it here before, moves with every node and every new triangulation, and
#: the mesh never settled. Over 1,522 offshore and other farms 0.95 settles
#: all but one within 250 iterations, that one in 317, and keeps the closest
#: pair at 0.98 of the median spacing on average. 0.9 settles sooner but
#: packs the nodes 3 % closer, and 1.0 and more press the mesh against its
#: edge, where interior nodes wedging between boundary nodes creep on for
#: hundreds of iterations.
_REST_LENGTH = 0.95
#: Most of its last move a node carries into the next one, which speeds the
#: slowest moves, a group of nodes sliding together. The share grows as
#: ``k / (k + 3)`` over the ``k`` iterations since the forces last turned
#: against the moves, when it starts again from none, the gradient restart
#: of O'Donoghue and Candès (2015). A fixed 0.7 settled 600 offshore and
#: other farms in 65 iterations on average and 248 for 99 in 100. Growing to
#: 0.9 and restarting settles them in 53 and 156.
_MOMENTUM = 0.9
#: Largest share of the area cutting off one sharp tip may remove.
_MAX_TIP_SHARE = 0.05
#: Share of the target spacing a node may move after a triangulation before
#: the nodes are triangulated again, the ``ttol`` of DistMesh. Between
#: triangulations the edges stay fixed while the nodes keep moving. Over
#: the 600 farms 0.2 keeps the closest pairs and settling of 0.1, and
#: triangulates a farm of 1,341 turbines 9 times, not 16.
_RETRIANGULATION = 0.2
def _sample_points(geometry, count, rng):
minx, miny, maxx, maxy = geometry.bounds
points = np.empty((0, 2))
attempts = 0
maximum_attempts = max(10_000, 2_000 * count)
while len(points) < count and attempts < maximum_attempts:
batch_size = min(max(4 * (count - len(points)), 64), 4096)
candidates = np.column_stack(
[
rng.uniform(minx, maxx, batch_size),
rng.uniform(miny, maxy, batch_size),
]
)
covered = candidates[shapely.covers(geometry, shapely.points(candidates))]
points = np.concatenate([points, covered[: count - len(points)]])
attempts += batch_size
if len(points) != count:
raise ValueError("could not initialize points inside the boundary")
return points
def _boundary_rings(geometry):
return (geometry.exterior, *geometry.interiors)
def _boundary_points(geometry, count, rng):
if count == 0:
return np.empty((0, 2)), []
rings = _boundary_rings(geometry)
counts = _allocate_counts(count, [ring.length for ring in rings])
points = []
ring_counts = []
for ring, ring_count in zip(rings, counts):
ring_count = int(ring_count)
ring_counts.append(ring_count)
if ring_count:
points.append(
_contained_perimeter_points(
geometry,
ring_count,
float(rng.random()),
line=ring,
)
)
combined = np.vstack(points) if points else np.empty((0, 2))
return combined, ring_counts
def _truncated_tips(geometry, width):
"""``geometry`` with each convex corner under 60 degrees cut where it is ``width`` wide.
Nodes either side of a sharp tip are closer than the spacing wherever
the wedge is narrower than it, so the mesh starts where the wedge is one
spacing wide, between two corners of at least 90 degrees. Returns the
shape and the number of corners cut. A cut that would remove more than
a twentieth of the area, or leave the shape invalid or split, is not made.
"""
coordinates = np.asarray(geometry.exterior.coords)[:-1]
cuts = []
for index, tip in enumerate(coordinates):
before = coordinates[index - 1] - tip
after = coordinates[(index + 1) % len(coordinates)] - tip
lengths = np.linalg.norm(before), np.linalg.norm(after)
if min(lengths) == 0.0:
continue
before, after = before / lengths[0], after / lengths[1]
angle = float(np.arccos(np.clip(before @ after, -1.0, 1.0)))
if angle >= np.pi / 3.0:
continue
bisector = (before + after) / np.linalg.norm(before + after)
probe = tip + 1e-3 * min(lengths) * bisector
# A reflex corner is as sharp seen from outside, and is left alone.
if not geometry.covers(shapely.Point(*probe)):
continue
# A rectangle from the cut line out past the tip and both edges,
# which it crosses at 60 degrees or more. A triangle on the tip and
# the edges themselves left, by rounding, a sliver along an edge, or
# the cut as a hole joined to the tip by one.
across = np.array([-bisector[1], bisector[0]])
depth = 0.5 * width / np.tan(0.5 * angle)
margin = 0.01 * width
cut = shapely.Polygon(
[
tip + depth * bisector + (0.5 * width + margin) * across,
tip + depth * bisector - (0.5 * width + margin) * across,
tip - margin * bisector - (0.5 * width + margin) * across,
tip - margin * bisector + (0.5 * width + margin) * across,
]
)
if cut.intersection(geometry).area <= _MAX_TIP_SHARE * geometry.area:
cuts.append(cut)
if not cuts:
return geometry, 0
truncated = geometry.difference(shapely.union_all(cuts))
if truncated.geom_type != "Polygon" or not truncated.is_valid:
return geometry, 0
return truncated, len(cuts)
def _initial_interior_points(geometry, count, boundary_points, rng):
if count == 0:
return np.empty((0, 2))
candidates = _sample_points(geometry, max(4 * count, count + 32), rng)
return _spread_selection(candidates, count, boundary_points, rng)
def _ring_segments(ring):
"""The starts, unit directions and arc-length ends of ``ring``'s segments.
The arc lengths run from 0 at the first vertex to the ring's length.
"""
coordinates = np.asarray(ring.coords)
sides = np.diff(coordinates, axis=0)
lengths = np.hypot(sides[:, 0], sides[:, 1])
ends = np.cumsum(lengths)
# A repeated vertex leaves a segment no arc length can land inside.
kept = lengths > 0
return (
coordinates[:-1][kept],
sides[kept] / lengths[kept, None],
np.concatenate([[0.0], ends[:-1]])[kept],
ends[kept],
)
def _ring_points(along, segments):
"""The points at arc lengths ``along`` of a ring."""
starts, directions, arcs, _ = segments
index = np.searchsorted(arcs, along, side="right") - 1
return starts[index] + (along - arcs[index])[:, None] * directions[index]
def _slid(along, moves, segments):
"""Arc lengths ``along`` advanced by the part of ``moves`` along the ring.
A node at a vertex moves on along whichever segment ``moves`` leads it
into. It stops at a vertex where the segment ahead would lead it back, as
either side of a corner, where it would otherwise cross the corner to
and fro and never settle.
"""
_, directions, arcs, ends = segments
index = np.searchsorted(arcs, along, side="right") - 1
ahead = np.sum(moves * directions[index], axis=1)
behind = np.sum(moves * directions[index - 1], axis=1)
after = np.sum(moves * directions[(index + 1) % len(arcs)], axis=1)
at_vertex = along == arcs[index]
slide = np.where(at_vertex & (ahead < 0.0), np.minimum(behind, 0.0), ahead)
slid = along + slide
slid = np.where((slid > ends[index]) & (after < 0.0), ends[index], slid)
slid = np.where(
~at_vertex & (slid < arcs[index]) & (behind > 0.0), arcs[index], slid
)
return slid % ends[-1]
def _triangulation(positions):
"""The Delaunay triangles of ``positions``."""
from scipy.spatial import Delaunay as _ScipyDelaunay
from scipy.spatial import QhullError
try:
return _ScipyDelaunay(positions).simplices
except QhullError as exc:
raise ValueError("could not triangulate the turbine positions") from exc
def _edges(triangles, count):
"""The sides of ``triangles`` over ``count`` nodes, once each and sorted."""
edges = np.vstack(
[
triangles[:, [0, 1]],
triangles[:, [1, 2]],
triangles[:, [2, 0]],
]
)
# One integer per sorted edge. Sorting the keys gives the row order of
# numpy.unique(axis=0) far more cheaply.
keys = np.unique(edges.min(axis=1) * count + edges.max(axis=1))
return np.column_stack([keys // count, keys % count])
def _mesh(positions, geometry):
"""The Delaunay triangles of ``positions`` that ``geometry`` covers, and their edges."""
if len(positions) < 3:
return np.empty((0, 3), dtype=int), np.empty((0, 2), dtype=int)
simplices = _triangulation(positions)
corners = positions[simplices]
sides = corners[:, 1:] - corners[:, :1]
flat = sides[:, 0, 0] * sides[:, 1, 1] - sides[:, 0, 1] * sides[:, 1, 0] == 0.0
keep = ~flat & shapely.covers(geometry, shapely.polygons(corners))
triangles = np.asarray(simplices[keep], dtype=int).reshape(-1, 3)
if not len(triangles):
return triangles, np.empty((0, 2), dtype=int)
return triangles, _edges(triangles, len(positions))
def _keep_distinct(previous, updated):
"""``updated``, with moves that land on another node undone.
Projection maps every node pushed past a corner onto the same point, and
the triangulation then drops the duplicates, so they would never separate.
``previous`` holds distinct nodes, so undoing moves always terminates.
"""
# Complex numbers sort by x and then y, and equal ones are equal nodes, so
# one sort shows most iterations have none, three times faster than the
# two keys of the full check.
ordered = np.sort(updated[:, 0] + 1j * updated[:, 1])
if not np.any(ordered[1:] == ordered[:-1]):
return updated
while True:
order = np.lexsort((updated[:, 1], updated[:, 0]))
ordered = updated[order]
same = np.all(ordered[1:] == ordered[:-1], axis=1)
if not np.any(same):
return updated
shared = np.zeros(len(updated), dtype=bool)
shared[order[1:][same]] = True
shared[order[:-1][same]] = True
undo = shared & np.any(updated != previous, axis=1)
if not np.any(undo):
return updated
updated[undo] = previous[undo]
[docs]
class Delaunay(LayoutAlgorithm):
"""Generate a boundary-conforming unstructured triangular layout.
Boundary nodes start at equal arc-length intervals and slide along their
ring. All nodes are connected by a Delaunay triangulation and repeatedly
pushed apart by edge forces, adapted from DistMesh (Persson and Strang
2004). ``boundary_count=None`` derives the number of boundary nodes from
the equilateral spacing of ``n_turbines`` triangles over the area, and a
node closer than 0.8 of the interval to an earlier one, across a sharp
corner or a narrow neck, becomes an interior node instead. Unlike
DistMesh, an edge repels its nodes while it is shorter than 0.95 of that
spacing, a fixed length rather than one scaled with the edges, and every
edge of the triangulation counts, also those leaving the shape. Each
iteration moves a node by ``relaxation`` times its mean edge force plus
up to 0.9 of its last move, a share that grows with the iterations since
the forces last turned against the moves and starts again from none
then. Iteration stops once no node moves more than ``tolerance`` times
that spacing, over 600 offshore and other farms within 156 iterations for
99 in 100 and 188 at most, so ``max_iterations`` defaults to 400 with room
to spare. The nodes are triangulated again once one has moved more
than 0.2 of that spacing since the last triangulation. On a multipart
boundary the parts share an explicit ``boundary_count`` in proportion to
their turbine counts.
"""
def __init__(
self,
n_turbines: int | None = None,
*,
density: dict[str, float] | None = None,
boundary_count: int | None = None,
max_iterations: int = 400,
tolerance: float = 1e-3,
relaxation: float = 1.0,
seed: int | np.random.Generator | None = None,
):
if boundary_count is not None:
boundary_count = _validate_integer(
"boundary_count", boundary_count, minimum=0
)
self.boundary_count = boundary_count
self.max_iterations = _validate_integer("max_iterations", max_iterations)
self.tolerance = _validate_positive("tolerance", tolerance)
self.relaxation = _validate_positive("relaxation", relaxation)
if self.relaxation > 1:
raise ValueError("relaxation must be in (0, 1]")
super().__init__(
n_turbines,
density=density,
seed=seed,
)
def _resolved_boundary_count(self, geometry, target_spacing):
n_turbines = self._required_n
if self.boundary_count is not None:
if self.boundary_count > n_turbines:
raise ValueError("boundary_count cannot exceed n_turbines")
return self.boundary_count
perimeter = sum(ring.length for ring in _boundary_rings(geometry))
estimate = int(round(perimeter / target_spacing))
minimum = min(3, n_turbines)
return min(n_turbines, max(minimum, estimate))
def _allocate_parts(self, areas):
counts, records = super()._allocate_parts(areas)
if self.boundary_count is None:
return counts, records
if self.boundary_count > self._required_n:
raise ValueError("boundary_count cannot exceed n_turbines")
# Shares in proportion to the counts never exceed a part's count.
shares = _allocate_counts(self.boundary_count, counts)
return counts, [{"boundary_count": int(share)} for share in shares]
def _part_layout(self, count, record):
child = super()._part_layout(count, record)
if "boundary_count" in record:
child.boundary_count = record["boundary_count"]
return child
def _fill_estimate(
self, boundary, spacing: float, mean: bool = False
) -> int | None:
"""A first count from one mesh at most of the turbines that fit.
A settled mesh keeps its closest pair near ``_CLOSEST_SHARE`` of its
median nearest-neighbor distance, so ``n (share median / spacing)^2``
from one draw near the answer lands within a few counts. The closest
pair itself hangs on one pair, which a narrow neck can pull a quarter
lower, and the search then stops among the counts that fit below. A
mean spacing sits near the median, so it takes no share.
"""
hexagonal = 2.0 * float(boundary.area) / (np.sqrt(3.0) * spacing**2)
trial_n = max(16, int(_TRIAL_SHARE * hexagonal))
try:
trial = self._with_count(trial_n).generate(boundary)
except ValueError:
return None
median = float(np.median(_nearest_neighbor_distances(trial.positions)))
if not np.isfinite(median):
return None
share = 1.0 if mean else _CLOSEST_SHARE
return int(trial_n * (share * median / spacing) ** 2)
def _generate_single(self, boundary, rng, seed):
target_spacing = float(
np.sqrt(2.0 * boundary.area / (np.sqrt(3.0) * self.n_turbines))
)
geometry, truncated = _truncated_tips(boundary.geometry, target_spacing)
# Each iteration tests every interior node against the boundary.
shapely.prepare(geometry)
boundary_count = self._resolved_boundary_count(geometry, target_spacing)
boundary_points, ring_counts = _boundary_points(geometry, boundary_count, rng)
if self.boundary_count is None and boundary_count > 1:
perimeter = sum(ring.length for ring in _boundary_rings(geometry))
kept = _kept_apart_mask(
boundary_points, _CROWDED_BOUNDARY * perimeter / boundary_count
)
boundary_points = boundary_points[kept]
ring_counts = [
int(part.sum()) for part in np.split(kept, np.cumsum(ring_counts)[:-1])
]
boundary_count = len(boundary_points)
interior_count = self._required_n - boundary_count
interior_points = _initial_interior_points(
geometry, interior_count, boundary_points, rng
)
positions = np.vstack([boundary_points, interior_points])
interior_geometry = _inset_interior(geometry)
# Each boundary node keeps its arc length along its own ring.
rings = []
along = np.empty(boundary_count)
start = 0
for ring, count in zip(_boundary_rings(geometry), ring_counts):
if count:
nodes = slice(start, start + count)
segments = _ring_segments(ring)
along[nodes] = (
shapely.line_locate_point(ring, shapely.points(positions[nodes]))
% segments[3][-1]
)
positions[nodes] = _ring_points(along[nodes], segments)
rings.append((nodes, segments))
start += count
rest_length = _REST_LENGTH * target_spacing
retriangulation = (_RETRIANGULATION * target_spacing) ** 2
converged = interior_count == 0
maximum_displacement = 0.0
iteration = 0
velocity = np.zeros_like(positions)
# Positions at the last triangulation, whose edges are kept until
# some node has moved far enough from them.
anchor = None
# Iterations since the forces last turned against the moves.
streak = 0
for iteration in range(1, self.max_iterations + 1):
if anchor is None or (
np.sum((positions - anchor) ** 2, axis=1).max() > retriangulation
):
if interior_count == 0 or len(positions) < 3:
break
# Every edge counts, those leaving the shape too, so nodes
# either side of a narrow notch push each other apart.
edges = _edges(_triangulation(positions), len(positions))
anchor = positions
ends = edges.T.ravel()
degrees = None
vectors = positions[edges[:, 0]] - positions[edges[:, 1]]
lengths = np.hypot(vectors[:, 0], vectors[:, 1])
usable = lengths > 0
if not np.all(usable):
if not np.any(usable):
break
edges = edges[usable]
vectors = vectors[usable]
lengths = lengths[usable]
ends = edges.T.ravel()
degrees = None
if degrees is None:
degrees = np.bincount(ends, minlength=len(positions)).astype(float)
degrees[degrees == 0] = 1.0
# Edges shorter than the rest length repel, so the truss only
# pushes nodes apart.
magnitudes = np.maximum(rest_length - lengths, 0.0)
edge_forces = vectors / lengths[:, None] * magnitudes[:, None]
# Each edge pushes its first node along its force and its second
# against it, summed in edge order.
signed = np.concatenate([edge_forces, -edge_forces])
forces = np.column_stack(
[
np.bincount(ends, weights=signed[:, 0], minlength=len(positions)),
np.bincount(ends, weights=signed[:, 1], minlength=len(positions)),
]
)
if float(np.sum(forces * velocity)) < 0.0:
streak = 0
momentum = min(_MOMENTUM, streak / (streak + 3.0))
streak += 1
moves = self.relaxation * forces / degrees[:, None] + momentum * velocity
updated = positions + moves
# Boundary nodes slide along their own ring, so the interior
# pressing on them spreads them out rather than wedging between.
slid = along.copy()
for nodes, segments in rings:
slid[nodes] = _slid(along[nodes], moves[nodes], segments)
updated[nodes] = _ring_points(slid[nodes], segments)
updated[boundary_count:] = _project_inside(
geometry, interior_geometry, updated[boundary_count:]
)
updated = _keep_distinct(positions, updated)
# A boundary node whose move was undone keeps its arc length.
undone = np.all(
updated[:boundary_count] == positions[:boundary_count], axis=1
)
along = np.where(undone, along, slid)
velocity = updated - positions
maximum_displacement = float(
np.sqrt(np.sum(velocity**2, axis=1).max(initial=0.0))
)
positions = updated
if maximum_displacement <= self.tolerance * target_spacing:
converged = True
break
# A cut corner re-nodes the shape, whose shared edges can then lie a
# rounding error outside the boundary, so the nodes are kept in it.
outline = boundary.geometry
positions = _project_inside(outline, _inset_interior(outline), positions)
triangles, edges = _mesh(positions, outline)
fit = {
"method": "delaunay",
"mesh_method": "distmesh",
"boundary_count": boundary_count,
"interior_count": interior_count,
"ring_counts": ring_counts,
"truncated_corners": truncated,
"target_spacing": target_spacing,
"iterations": iteration,
"converged": converged,
"maximum_displacement": maximum_displacement,
"tolerance": self.tolerance,
"relaxation": self.relaxation,
"triangles": triangles.tolist(),
"edges": edges.tolist(),
}
return self._single_result(positions, boundary, seed, fit)
