"""Modifiers of boundaries, string guides and generated layouts.
Layout modifiers change a finished farm and attach to a layout with
:meth:`~oowl.layouts.LayoutAlgorithm.modify`. Boundary and string modifiers
are re-exported so every modifier is available from one module.
"""
from __future__ import annotations
import math
import typing
from abc import ABC, abstractmethod
from collections.abc import Callable
from copy import copy
from functools import lru_cache
from itertools import chain
import numpy as np
import shapely
from ._random import pin_seed
from ._repr import _ArgumentRepr
from ._validation import (
_validate_integer,
_validate_nonnegative,
_validate_positive,
_validate_real,
)
from .boundaries.modifiers import BoundaryModifier, ExclusionZones, Setback
from .boundaries.partition import Partition
from .layouts.base import Farm, _median_spacing
from .strings.modifiers import Breaks, Meander, StringModifier
__all__ = [
"LayoutModifier",
"BoundaryModifier",
"StringModifier",
"ExclusionZones",
"Setback",
"Partition",
"Breaks",
"Meander",
"PositionNoise",
"Relaxation",
"TurbineRemoval",
]
[docs]
class LayoutModifier(_ArgumentRepr, ABC):
"""Transform a completed layout after its base algorithm has run."""
seed = None
def _placed_for(self, n_kept: int) -> int:
"""Turbines to place for ``n_kept`` to remain after this modifier.
A layout given ``n_turbines`` places them so that many remain.
"""
return n_kept
def _held_apart(self, spacing: float) -> "LayoutModifier":
"""This modifier, bringing no pair of turbines closer than ``spacing``.
A layout filled at a spacing runs its modifiers once on the filled
farm through this. A modifier that only removes turbines returns
itself, and one that moves them returns a copy held to ``spacing``.
A layout raises ``ValueError`` when a pair ends up closer anyway.
"""
return self
[docs]
@abstractmethod
def apply(self, layout: Farm, rng: np.random.Generator, seed: int) -> Farm:
"""Return a new ``Farm``, leaving ``layout`` untouched.
``rng`` is seeded from ``seed``, this modifier's effective seed.
Implementations append a record to ``metadata["modifiers"]``.
"""
def _own_pieces(layout: Farm) -> list | None:
"""The boundary polygon each turbine starts in, or ``None`` without one."""
if layout.boundary is None:
return None
geometry = layout.boundary.geometry
polygons = list(getattr(geometry, "geoms", [geometry]))
if len(polygons) == 1:
return polygons * len(layout.positions)
tree = shapely.STRtree(polygons)
points = shapely.points(layout.positions)
owners = [polygons[0]] * len(points)
for point, polygon in zip(*tree.query(points, predicate="intersects")):
owners[int(point)] = polygons[int(polygon)]
return owners
def _updated_ranges(metadata, kept):
"""Update contiguous multipart slices after subsetting positions."""
for key in ("components", "parts"):
if key not in metadata:
continue
records = []
new_start = 0
for original in metadata[key]:
record = dict(original)
old_start = int(record["start"])
old_stop = int(record["stop"])
count = int(np.count_nonzero((kept >= old_start) & (kept < old_stop)))
record["n_turbines"] = count
record["start"] = new_start
record["stop"] = new_start + count
records.append(record)
new_start += count
metadata[key] = records
if key == "parts" and "counts" in metadata:
metadata["counts"] = [record["n_turbines"] for record in records]
def _subset_layout(layout, kept, modifier_record):
kept = np.sort(np.asarray(kept, dtype=int))
metadata = dict(layout.metadata)
_updated_ranges(metadata, kept)
if "origin" in metadata:
metadata["origin"] = np.asarray(metadata["origin"])[kept].tolist()
metadata["achieved_n"] = int(len(kept))
modifiers = list(metadata.get("modifiers", []))
modifiers.append(modifier_record)
metadata["modifiers"] = modifiers
turbines = layout.turbines
if isinstance(turbines, tuple):
turbines = tuple(turbines[index] for index in kept)
return Farm(
layout.positions[kept],
layout.boundary,
metadata=metadata,
guide=layout.guide,
resource=layout.resource,
turbines=turbines,
surface=layout.surface,
)
def _grid_rows(positions, basis, row_axis):
"""Return lattice-ordered position indices and row-start indicators."""
positions = np.asarray(positions, dtype=float)
basis = np.asarray(basis, dtype=float)
if basis.shape != (2, 2) or not np.all(np.isfinite(basis)):
raise ValueError("grid metadata contains an invalid grid_basis")
try:
coordinates = (positions - positions[0]) @ np.linalg.inv(basis)
except np.linalg.LinAlgError as exc:
raise ValueError("grid metadata contains a singular grid_basis") from exc
lattice = np.rint(coordinates).astype(int)
if not np.allclose(coordinates, lattice, atol=1e-5, rtol=0.0):
raise ValueError("positions do not align with their recorded grid_basis")
other_axis = 1 - row_axis
order = np.lexsort((lattice[:, row_axis], lattice[:, other_axis]))
ordered_rows = lattice[order, other_axis]
row_starts = np.r_[True, ordered_rows[1:] != ordered_rows[:-1]]
return order, row_starts
def _markov_order(layout, row_axis):
"""Find row order for a single or multipart Grid result."""
return _recorded_rows(layout.metadata, layout.positions, row_axis)
def _recorded_rows(metadata, positions, row_axis):
"""Row order of ``positions`` laid out as ``metadata`` records.
A multipart record holds a component per part, and a partition cell in
several polygons is a multipart record of its own, so components nest.
Each component's ``start`` and ``stop`` index its parent's positions.
"""
if metadata.get("kind") == "single":
basis = metadata.get("fit", {}).get("grid_basis")
if basis is None:
raise ValueError("markov turbine removal requires a Grid")
return _grid_rows(positions, basis, row_axis)
if metadata.get("kind") == "multipart":
orders = []
starts = []
for component in metadata.get("components", []):
start, stop = int(component["start"]), int(component["stop"])
if stop <= start:
continue
order, row_starts = _recorded_rows(
component.get("layout") or {}, positions[start:stop], row_axis
)
orders.append(order + start)
starts.append(row_starts)
if orders:
return np.concatenate(orders), np.concatenate(starts)
raise ValueError("markov turbine removal requires a Grid")
def _sample_markov_mask(n_positions, n_remove, row_starts, penalty, rng):
"""Sample a removal mask with exactly ``n_remove`` removed positions.
Each mask has weight ``exp(-penalty * transitions)``, where transitions
count kept and removed neighbors within a row. Forward filtering then
backward sampling draws a mask exactly from that distribution.
"""
row_starts = np.asarray(row_starts, dtype=bool)
alpha = np.full((n_positions, n_remove + 1, 2), -np.inf)
alpha[0, 0, 0] = 0.0
if n_remove:
alpha[0, 1, 1] = 0.0
for index in range(1, n_positions):
previous = alpha[index - 1]
if row_starts[index]:
alpha[index, :, 0] = np.logaddexp(previous[:, 0], previous[:, 1])
alpha[index, 1:, 1] = np.logaddexp(previous[:-1, 0], previous[:-1, 1])
else:
alpha[index, :, 0] = np.logaddexp(previous[:, 0], previous[:, 1] - penalty)
alpha[index, 1:, 1] = np.logaddexp(
previous[:-1, 0] - penalty, previous[:-1, 1]
)
mask = np.zeros(n_positions, dtype=bool)
remaining = n_remove
logits = alpha[-1, remaining]
probabilities = np.exp(logits - np.max(logits))
probabilities /= probabilities.sum()
state = int(rng.choice(2, p=probabilities))
mask[-1] = bool(state)
for index in range(n_positions - 1, 0, -1):
previous_removed = remaining - state
logits = alpha[index - 1, previous_removed].copy()
if not row_starts[index]:
logits -= penalty * (np.arange(2) != state)
probabilities = np.exp(logits - np.max(logits))
probabilities /= probabilities.sum()
state = int(rng.choice(2, p=probabilities))
mask[index - 1] = bool(state)
remaining = previous_removed
return mask
[docs]
class TurbineRemoval(LayoutModifier):
"""Remove turbines randomly or in correlated grid rows.
``n_remove`` removes that many turbines, and a layout given
``n_turbines`` places them on top, so exactly ``n_turbines`` remain.
``relative_n_remove`` in ``[0, 1)`` removes that share of the turbines
the layout places, rounded down, so at least one always remains however
small a sampled farm turns out. A layout given ``n_turbines`` then places
the fewest turbines that leave exactly that many. Give one of the two.
``method="random"`` removes a uniform random subset. ``method="markov"``
needs a :class:`~oowl.layouts.Grid` layout and removes runs of neighbors
along its rows, which follow the first grid vector for ``row_axis=0`` and
the second for ``row_axis=1``. A larger ``transition_penalty`` gives
longer runs.
"""
def __init__(
self,
n_remove: int | None = None,
*,
relative_n_remove: float | None = None,
method: typing.Literal["random", "markov"] = "random",
transition_penalty: float | None = None,
row_axis: int = 0,
seed: int | np.random.Generator | None = None,
) -> None:
if (n_remove is None) == (relative_n_remove is None):
raise TypeError("give one of n_remove and relative_n_remove")
if n_remove is not None:
n_remove = _validate_integer("n_remove", n_remove, minimum=0)
if relative_n_remove is not None:
relative_n_remove = _validate_real("relative_n_remove", relative_n_remove)
if not 0.0 <= relative_n_remove < 1.0:
raise ValueError("relative_n_remove must be in [0, 1)")
self.n_remove = n_remove
self.relative_n_remove = relative_n_remove
if method not in ("random", "markov"):
raise ValueError("method must be 'random' or 'markov'")
if method == "markov" and transition_penalty is None:
raise ValueError("transition_penalty is required for method='markov'")
if transition_penalty is not None:
transition_penalty = _validate_nonnegative(
"transition_penalty", transition_penalty
)
if isinstance(row_axis, (bool, np.bool_)) or not isinstance(
row_axis, (int, float, np.integer, np.floating)
):
raise TypeError(f"row_axis must be 0 or 1, not {type(row_axis).__name__}")
if row_axis not in (0, 1):
raise ValueError("row_axis must be 0 or 1")
self.method = method
self.transition_penalty = transition_penalty
self.row_axis = int(row_axis)
self.seed = pin_seed(seed)
def _removed_from(self, n_placed: int) -> int:
"""Turbines this removes from a layout of ``n_placed``."""
share = self.relative_n_remove
if share is None:
return self.n_remove or 0
return math.floor(share * n_placed)
def _placed_for(self, n_kept: int) -> int:
share = self.relative_n_remove
if share is None:
return n_kept + (self.n_remove or 0)
# What remains grows by none or one turbine per turbine placed, so
# the fewest placed is found from the estimate by single steps.
placed = max(n_kept, math.ceil(n_kept / (1.0 - share)))
while placed - self._removed_from(placed) < n_kept:
placed += 1
while placed > n_kept and placed - 1 - self._removed_from(placed - 1) >= n_kept:
placed -= 1
return placed
[docs]
def apply(self, layout: Farm, rng: np.random.Generator, seed: int) -> Farm:
"""Return a new farm with the turbines this removes dropped.
Removal is uniform for ``method="random"`` and row-correlated for
``method="markov"``. Records the seed, method and turbine counts in
``metadata["modifiers"]`` and updates ``metadata["achieved_n"]``.
Raises ``ValueError`` if ``n_remove`` would leave no turbine.
"""
n_remove = self._removed_from(layout.n_turbines)
if n_remove >= layout.n_turbines:
raise ValueError(
f"n_remove={n_remove} must leave at least one of the "
f"{layout.n_turbines} turbines the layout placed, so remove "
"a share of them with relative_n_remove"
)
keep_count = layout.n_turbines - n_remove
record = {
"type": "turbine_removal",
"method": self.method,
"seed": int(seed),
"input_n": layout.n_turbines,
"removed_n": n_remove,
"output_n": keep_count,
}
if self.relative_n_remove is not None:
record["relative_n_remove"] = self.relative_n_remove
if self.method == "random":
kept = rng.choice(layout.n_turbines, size=keep_count, replace=False)
return _subset_layout(layout, kept, record)
order, row_starts = _markov_order(layout, self.row_axis)
removed = _sample_markov_mask(
layout.n_turbines,
n_remove,
row_starts,
self.transition_penalty,
rng,
)
kept = order[~removed]
transitions = int(
np.count_nonzero((removed[1:] != removed[:-1]) & ~row_starts[1:])
)
record.update(
transition_penalty=self.transition_penalty,
row_axis=self.row_axis,
row_transitions=transitions,
)
return _subset_layout(layout, kept, record)
[docs]
class PositionNoise(LayoutModifier):
"""Add boundary-constrained Gaussian turbine-position noise.
``relative_standard_deviation`` is a fraction of the layout's median
nearest-neighbor spacing, not a distance. When the layout carries a
:class:`~oowl.strings.StringGuide`, a turbine on it moves only along it,
so it stays on the string. Turbines off the guide, such as those a
:class:`~oowl.layouts.Mixture` takes from another source, move freely.
``minimum_spacing`` in meters rejects moves that come closer to another
turbine, and each turbine gets ``max_attempts`` draws.
"""
def __init__(
self,
relative_standard_deviation: float,
*,
minimum_spacing: float | None = None,
seed: int | np.random.Generator | None = None,
max_attempts: int = 100,
) -> None:
self.relative_standard_deviation = _validate_nonnegative(
"relative_standard_deviation", relative_standard_deviation
)
self.minimum_spacing = (
None
if minimum_spacing is None
else _validate_nonnegative("minimum_spacing", minimum_spacing)
)
self.seed = pin_seed(seed)
self.max_attempts = _validate_integer("max_attempts", max_attempts)
def _held_apart(self, spacing: float) -> "PositionNoise":
held = copy(self)
held.minimum_spacing = max(self.minimum_spacing or 0.0, spacing)
return held
@staticmethod
def _reference_spacing(layout: Farm) -> float:
"""Median nearest-neighbor spacing, else the spacing of the area."""
median = _median_spacing(layout.positions)
if median is not None:
return median
if layout.boundary is not None:
return float(np.sqrt(layout.boundary.area / max(layout.n_turbines, 1)))
return 1.0
[docs]
def apply(self, layout: Farm, rng: np.random.Generator, seed: int) -> Farm:
"""Return a new farm with Gaussian noise added to every position.
The noise has ``relative_standard_deviation`` times the median
nearest-neighbor spacing. Turbines move in order, each trying up to
``max_attempts`` draws against the current positions of the others.
A turbine with no boundary- and spacing-feasible candidate stays put.
A turbine never leaves the piece of the boundary it starts in. Raises
``ValueError`` when the input already violates ``minimum_spacing``.
"""
positions = np.asarray(layout.positions, dtype=float).copy()
# The margin covers rounding in the distances a fill measured.
if self.minimum_spacing is not None and layout.minimum_spacing() < (
self.minimum_spacing * (1.0 - 1e-9)
):
raise ValueError("layout does not satisfy the requested minimum_spacing")
reference_spacing = self._reference_spacing(layout)
sigma = self.relative_standard_deviation * reference_spacing
moved = 0
along_string = False
pieces = _own_pieces(layout)
if sigma > 0:
from .strings import StringGuide
n_turbines = len(positions)
on_string = np.zeros(n_turbines, dtype=bool)
guide = None
nearest = np.zeros(n_turbines, dtype=int)
distances = np.zeros(n_turbines)
if isinstance(layout.guide, StringGuide):
guide = _GuideLines(layout.guide.geometry)
nearest, distances, offsets = guide.project(positions)
# A mixture's turbines from other sources lie off the guide
# and take ordinary noise. Turbines on it differ from it by
# rounding only, at most a few coordinate units in the last
# place.
on_string = offsets <= max(
1e-9 * reference_spacing,
64.0 * float(np.spacing(np.abs(positions).max(initial=0.0))),
)
along_string = bool(on_string.any())
strings = (on_string, guide, nearest, distances)
moved = _jitter(
positions,
_NormalDraws(rng, sigma, 2 * n_turbines),
self.max_attempts,
pieces,
strings,
None
if self.minimum_spacing is None or n_turbines == 1
else _MovingNeighbors(positions, self.minimum_spacing),
)
metadata = dict(layout.metadata)
modifiers = list(metadata.get("modifiers", []))
modifiers.append(
dict(
type="position_noise",
seed=int(seed),
relative_standard_deviation=self.relative_standard_deviation,
minimum_spacing=self.minimum_spacing,
reference_spacing=reference_spacing,
moved=moved,
along_string=along_string,
)
)
metadata["modifiers"] = modifiers
return Farm(
positions,
layout.boundary,
metadata=metadata,
guide=layout.guide,
resource=layout.resource,
turbines=layout.turbines,
surface=layout.surface,
)
class _GuideLines:
"""The lines of a string guide, measured once for many points along them."""
def __init__(self, geometry) -> None:
self.lines = shapely.get_parts(geometry)
self.lengths = shapely.length(self.lines)
self._vertices = []
self._segments = []
self._starts = []
for line in self.lines:
vertices = shapely.get_coordinates(line)
steps = np.diff(vertices, axis=0)
segments = np.sqrt(steps[:, 0] * steps[:, 0] + steps[:, 1] * steps[:, 1])
self._vertices.append(vertices)
self._segments.append(segments)
self._starts.append(np.concatenate(([0.0], np.cumsum(segments))))
def project(self, positions):
"""Each position's nearest line, distance along it and offset from it.
The first line wins a tie.
"""
points = shapely.points(positions)
nearest = np.zeros(len(points), dtype=int)
distances = np.zeros(len(points))
offsets = np.full(len(points), np.inf)
for index, line in enumerate(self.lines):
along = shapely.line_locate_point(line, points)
gaps = self.points(index, along) - positions
offset = np.sqrt(gaps[:, 0] * gaps[:, 0] + gaps[:, 1] * gaps[:, 1])
closer = offset < offsets
nearest[closer] = index
distances[closer] = along[closer]
offsets[closer] = offset[closer]
return nearest, distances, offsets
def points(self, index, distances):
"""The points ``distances`` along line ``index`` from its start.
They match ``shapely.line_interpolate_point`` bit for bit. A distance
at most zero gives the start, and one past the end gives the end.
"""
vertices = self._vertices[index]
segments = self._segments[index]
starts = self._starts[index]
last = len(segments) - 1
segment = np.minimum(np.searchsorted(starts[1:], distances, side="right"), last)
first = vertices[segment]
second = vertices[segment + 1]
fraction = ((distances - starts[segment]) / segments[segment])[:, np.newaxis]
points = (second - first) * fraction + first
points = np.where(
fraction <= 0.0, first, np.where(fraction >= 1.0, second, points)
)
points[distances >= starts[-1]] = vertices[-1]
points[distances <= 0.0] = vertices[0]
return points
class _NormalDraws:
"""The ``rng.normal(scale=sigma)`` draws in order, read ahead in blocks.
Separate ``rng.normal`` calls of any size continue one stream, so
reading ahead and consuming only the draws used gives the values of one
call per draw.
"""
def __init__(self, rng: np.random.Generator, sigma: float, block: int) -> None:
self._rng = rng
self._sigma = sigma
self._block = block
self._buffer = np.empty(0)
self._start = 0
def peek(self, count: int) -> np.ndarray:
"""The next ``count`` draws, leaving them unconsumed."""
if self._start + count > len(self._buffer):
self._buffer = np.concatenate(
[
self._buffer[self._start :],
self._rng.normal(scale=self._sigma, size=max(count, self._block)),
]
)
self._start = 0
return self._buffer[self._start : self._start + count]
def consume(self, count: int) -> None:
"""Mark the next ``count`` draws used."""
self._start += count
class _MovingNeighbors:
"""Spacing tests against ``positions`` while its turbines move one by one.
A tree holds the starting positions, and ``reach`` bounds how far any
turbine has moved since, so a search ``reach`` wider finds every turbine
near a candidate.
"""
def __init__(self, positions: np.ndarray, spacing: float) -> None:
from scipy.spatial import cKDTree # ty: ignore[unresolved-import]
self._positions = positions
self._start = positions.copy()
self._tree = cKDTree(self._start)
self._spacing = spacing
self._reach = 0.0
def _search(self, length, reach):
# The margin keeps every turbine whose distance rounds near the spacing.
return (length + self._spacing + reach) * (1.0 + 1e-9)
def clear_attempts(self, index: int, candidates: np.ndarray) -> np.ndarray:
"""Which ``candidates`` for unmoved turbine ``index`` keep the spacing."""
center = self._start[index]
offsets = candidates - center
length = float(np.sqrt(np.max((offsets * offsets).sum(axis=1))))
near = np.asarray(
self._tree.query_ball_point(center, self._search(length, self._reach)),
dtype=int,
)
near = near[near != index]
differences = self._positions[near] - candidates[:, np.newaxis]
squared = (differences * differences).sum(axis=2)
return np.sqrt(squared.min(axis=1, initial=np.inf)) >= self._spacing
def clear_first(self, turbines: np.ndarray, candidates: np.ndarray) -> np.ndarray:
"""Which ``candidates`` of consecutive unmoved ``turbines`` keep the spacing.
Each turbine sees the earlier ones of ``turbines`` at their candidates.
"""
offsets = candidates - self._start[turbines]
lengths = np.sqrt((offsets * offsets).sum(axis=1))
found = self._tree.query_ball_point(
self._start[turbines],
self._search(lengths, max(self._reach, float(lengths.max()))),
return_sorted=False,
)
counts = np.fromiter(map(len, found), dtype=int, count=len(found))
near = np.fromiter(
chain.from_iterable(found), dtype=int, count=int(counts.sum())
)
rows = np.repeat(np.arange(len(turbines)), counts)
kept = near != turbines[rows]
near, rows = near[kept], rows[kept]
others = self._positions[near]
earlier = (near >= turbines[0]) & (near < turbines[rows])
others[earlier] = candidates[near[earlier] - turbines[0]]
differences = others - candidates[rows]
closest = np.full(len(candidates), np.inf)
np.minimum.at(closest, rows, (differences * differences).sum(axis=1))
return np.sqrt(closest) >= self._spacing
def moved(self, turbines: np.ndarray) -> None:
"""Widen the search after ``turbines`` moved."""
if len(turbines):
offsets = self._positions[turbines] - self._start[turbines]
self._reach = max(
self._reach, float(np.sqrt((offsets * offsets).sum(axis=1).max()))
)
def _jitter(positions, draws, max_attempts, pieces, strings, neighbors) -> int:
"""Move each turbine in order to its first feasible noisy candidate.
``strings`` holds which turbines move along the guide, its
:class:`_GuideLines` or ``None``, and each turbine's nearest line and
distance along it. Returns how many turbines moved. A block of turbines
tries its first draws at once, each seeing the earlier ones at their
candidates, which holds up to the first that fails. That turbine tries
its remaining draws alone, and the next block starts after it. While
first draws keep failing, each turbine tries all its draws alone.
"""
on_string, guide, nearest, distances = strings
widths = np.where(on_string, 1, 2)
def along_guide(turbines, noise):
"""Points ``noise`` along the guide from on-string ``turbines``."""
points = np.empty((len(turbines), 2))
for line in np.unique(nearest[turbines]):
rows = nearest[turbines] == line
points[rows] = guide.points(
line,
np.clip(
distances[turbines[rows]] + noise[rows],
0.0,
guide.lengths[line],
),
)
return points
if pieces is not None:
pieces = np.asarray(pieces, dtype=object)
shapely.prepare(pieces)
def trial(turbines, noise, starts):
"""Candidates for ``turbines`` from ``noise`` at ``starts``, and which fit."""
candidates = np.empty((len(turbines), 2))
free = ~on_string[turbines]
candidates[free] = positions[turbines[free]] + np.column_stack(
(noise[starts[free]], noise[starts[free] + 1])
)
along = on_string[turbines]
if along.any():
candidates[along] = along_guide(turbines[along], noise[starts[along]])
if pieces is None:
return candidates, np.ones(len(turbines), dtype=bool)
return candidates, shapely.intersects_xy(
pieces[turbines], candidates[:, 0], candidates[:, 1]
)
def alone(index, tried):
"""Try turbine ``index`` from draw ``tried`` on, giving the one it took.
Returns -1 when no draw fits and the turbine stays put.
"""
width = int(widths[index])
count = 16
while tried < max_attempts:
count = min(count, max_attempts - tried)
noise = draws.peek(count * width)
if width == 1:
candidates = along_guide(np.full(count, index), noise)
else:
candidates = positions[index] + noise.reshape(count, 2)
feasible = (
np.ones(count, dtype=bool)
if pieces is None
else shapely.intersects_xy(
pieces[index], candidates[:, 0], candidates[:, 1]
)
)
if neighbors is not None and feasible.any():
feasible[feasible] = neighbors.clear_attempts(
index, candidates[feasible]
)
accepted = np.flatnonzero(feasible)
if len(accepted):
first = int(accepted[0])
draws.consume((first + 1) * width)
positions[index] = candidates[first]
if neighbors is not None:
neighbors.moved(np.array([index]))
return tried + first
draws.consume(count * width)
tried += count
count *= 4
return -1
moved = 0
start = 0
block = 4
while start < len(positions):
if block < 2:
# First draws mostly fail, so turbines go alone until one fits.
taken = alone(start, 0)
moved += taken >= 0
block = 4 if taken == 0 else 0
start += 1
continue
turbines = np.arange(start, min(start + block, len(positions)))
ends = np.cumsum(widths[turbines])
candidates, feasible = trial(
turbines, draws.peek(int(ends[-1])), ends - widths[turbines]
)
if neighbors is not None:
feasible &= neighbors.clear_first(turbines, candidates)
failed = np.flatnonzero(~feasible)
run = int(failed[0]) if len(failed) else len(turbines)
positions[turbines[:run]] = candidates[:run]
if neighbors is not None:
neighbors.moved(turbines[:run])
moved += run
start += run
block = min(256, 2 * run)
if run == len(turbines):
draws.consume(int(ends[-1]))
continue
# The failed turbine's first draw is used up too.
draws.consume(int(ends[run]))
moved += alone(start, 1) >= 0
start += 1
return moved
def _zero_gradient(positions):
"""Return a callable giving the zero position gradient."""
return lambda: np.zeros_like(positions)
def _pair_gradient(first, second, pair_gradient, n_turbines):
"""Accumulate equal and opposite pair contributions at each turbine."""
return np.column_stack(
[
np.bincount(first, weights=pair_gradient[:, axis], minlength=n_turbines)
- np.bincount(second, weights=pair_gradient[:, axis], minlength=n_turbines)
for axis in range(2)
]
)
# Pairs farther apart than this many spacings hold less than exp(-16) of the
# Gaussian energy, so leaving them out changes it by about 1e-7.
_GAUSSIAN_RADIUS = 4.0
def _gaussian_energy(positions, spacing, rng):
"""Return the Gaussian and spacing pair energy and a gradient callable.
The Gaussian term covers the pairs within ``_GAUSSIAN_RADIUS`` spacings,
which one KD-tree query finds, and the spacing penalty those closer than
``spacing``. ``rng`` gives exactly coincident pairs a direction.
"""
n_turbines = len(positions)
if n_turbines < 2:
return 0.0, _zero_gradient(positions)
from scipy.spatial import cKDTree # ty: ignore[unresolved-import]
pairs = cKDTree(positions).query_pairs(
_GAUSSIAN_RADIUS * spacing, output_type="ndarray"
)
first, second = pairs[:, 0], pairs[:, 1]
differences = positions[first] - positions[second]
squared_distance = np.einsum("ij,ij->i", differences, differences)
gaussian = np.exp(-squared_distance / spacing**2)
close = np.flatnonzero(squared_distance < spacing**2)
distance = np.sqrt(squared_distance[close])
violation = 1.0 - distance / spacing
energy = float(np.sum(gaussian)) / n_turbines + 0.5 * float(np.sum(violation**2))
# Randomize direction for coincident pairs.
directions = np.empty((len(close), 2))
coincident = distance <= spacing * np.finfo(float).eps
distinct = ~coincident
directions[distinct] = differences[close[distinct]] / distance[distinct, None]
if np.any(coincident):
angles = rng.uniform(0.0, 2.0 * np.pi, int(np.count_nonzero(coincident)))
directions[coincident] = np.column_stack((np.cos(angles), np.sin(angles)))
def gradient():
"""Return the Gaussian and spacing gradient at ``positions``."""
scale = (-2.0 / (n_turbines * spacing**2)) * gaussian
pair_gradient = scale[:, None] * differences
pair_gradient[close] += (-violation / spacing)[:, None] * directions
return _pair_gradient(first, second, pair_gradient, n_turbines)
return energy, gradient
@lru_cache(maxsize=8)
def _pair_indices(n_turbines):
"""Read-only upper-triangle index pairs, shared by every energy evaluation."""
first, second = np.triu_indices(n_turbines, k=1)
first.setflags(write=False)
second.setflags(write=False)
return first, second
# Pairs per block. Block arrays stay in cache, so the cost per pair does not
# grow with the turbine count.
_PAIR_BLOCK = 32768
def _pair_squared_distances(positions, first, second):
"""Yield slices of the pairs ``first``, ``second`` with their squared distances.
Each block holds at most ``_PAIR_BLOCK`` pairs. Its squared distances are
a fresh array that callers may overwrite.
"""
# A complex view gathers both coordinates of a turbine at once.
points = np.ascontiguousarray(positions, dtype=float).view(np.complex128).ravel()
for start in range(0, len(first), _PAIR_BLOCK):
block = slice(start, start + _PAIR_BLOCK)
differences = points[first[block]]
differences -= points[second[block]]
squared_distance = differences.real * differences.real
squared_distance += differences.imag * differences.imag
yield block, squared_distance
def _row_blocks(n_turbines):
"""Yield row ranges with about ``_PAIR_BLOCK`` upper-triangle pairs each.
Ranges hold at least two rows, except a final range with a single pair,
which keeps the reductions in ``_all_pair_gradient`` off the fast axis.
"""
start = 0
while start < n_turbines - 1:
rows = max(2, _PAIR_BLOCK // (n_turbines - start))
stop = min(n_turbines - 1, start + rows)
yield start, stop
start = stop
def _all_pair_gradient(positions, pair_scale, coincident=None):
r"""Return :math:`\sum_{j>i} w_{ij} - \sum_{k<i} w_{ki}` at each turbine ``i``.
The pair terms are :math:`w_{ij}=s_{ij}(p_i-p_j)`, where
``pair_scale(start, stop, differences)`` returns the scales of rows
``start`` to ``stop`` against columns ``start + 1`` onward. An optional
``coincident`` holds turbine index pairs and the gradient rows that
replace their terms.
Both sums run in the same order as ``np.bincount`` over upper-triangle
pairs, which keeps the result bit for bit. NumPy adds sequentially when
reducing along a slow axis, and here that axis vectorizes over turbines.
"""
n_turbines = len(positions)
x = np.ascontiguousarray(positions[:, 0], dtype=float)
y = np.ascontiguousarray(positions[:, 1], dtype=float)
first_sums = np.zeros((2, n_turbines))
second_sums = np.zeros((2, n_turbines))
for start, stop in _row_blocks(n_turbines):
rows, columns = stop - start, slice(start + 1, n_turbines)
buffer = np.empty((2, rows + 1, n_turbines - start - 1))
terms = buffer[:, 1:]
np.subtract(x[start:stop, None], x[None, columns], out=terms[0])
np.subtract(y[start:stop, None], y[None, columns], out=terms[1])
terms *= pair_scale(start, stop, terms)
if coincident is not None:
pairs, pair_gradient = coincident
inside = (pairs[:, 0] >= start) & (pairs[:, 0] < stop)
row, column = pairs[inside].T
terms[:, row - start, column - start - 1] = pair_gradient[inside].T
# Zero the columns at or before each row, which are not pairs.
terms[:, :, :rows][:, np.tri(rows, k=-1, dtype=bool)] = 0.0
# Each column continues its running sum over earlier rows.
buffer[:, 0] = second_sums[:, columns]
np.add.reduce(buffer, axis=1, out=second_sums[:, columns])
# Transposing puts each row's pairs along the slow axis. Starting
# from zero matches ``np.bincount`` for signed zeros.
np.add.reduce(
np.ascontiguousarray(terms.transpose(0, 2, 1)),
axis=1,
initial=0.0,
out=first_sums[:, start:stop],
)
return np.ascontiguousarray((first_sums - second_sums).T)
def _repulsion_energy(positions, spacing, rng, softening):
r"""Return dimensionless all-pairs inverse-square energy and a gradient callable.
``softening`` regularizes coincident points. The gradient is only computed
on request, since the line search never uses it for rejected candidates.
"""
n_turbines = len(positions)
if n_turbines < 2:
return 0.0, _zero_gradient(positions)
first, second = _pair_indices(n_turbines)
pair_energy = np.empty(len(first))
coincident_blocks = []
for block, denominator in _pair_squared_distances(positions, first, second):
coincident_blocks.append(np.flatnonzero(denominator == 0.0) + block.start)
denominator += softening**2
np.divide(spacing**2, denominator, out=pair_energy[block])
energy = float(np.sum(pair_energy))
# Randomize direction for exact coincidences.
coincident = np.concatenate(coincident_blocks)
if len(coincident):
angles = rng.uniform(0.0, 2.0 * np.pi, len(coincident))
def pair_scale(start, stop, differences):
"""Return the inverse-square gradient scales of a row block."""
scale = differences[0] * differences[0]
scale += differences[1] * differences[1]
scale += softening**2
np.square(scale, out=scale)
np.divide(-2.0 * spacing**2, scale, out=scale)
return scale
def gradient():
"""Return the inverse-square gradient at ``positions``."""
if not len(coincident):
return _all_pair_gradient(positions, pair_scale)
directions = np.column_stack((np.cos(angles), np.sin(angles)))
pairs = np.column_stack((first[coincident], second[coincident]))
pair_gradient = (-(spacing**2) / (2.0 * softening**3)) * directions
return _all_pair_gradient(positions, pair_scale, (pairs, pair_gradient))
return energy, gradient
def _limit_displacement(candidate, initial, maximum):
if maximum is None:
return candidate
displacement = candidate - initial
distance = np.linalg.norm(displacement, axis=1)
outside = distance > maximum
if np.any(outside):
displacement[outside] *= (maximum / distance[outside])[:, None]
return initial + displacement
# Line search, boundary bisection and the repulsion softening resolve
# distances to this fraction of the spacing.
_DISTANCE_TOLERANCE = 1e-3
# A relaxation stops once its energy falls by less than its tolerance over
# this many accepted steps.
_STALL_STEPS = 5
def _slid_to_edge(geometry, edge, target, distance_tolerance):
"""The edge point nearest each ``target``, ``distance_tolerance`` inside.
Rows are NaN for targets inside ``geometry`` or on its edge, which have
no inward direction to slide along.
"""
slid = np.full_like(target, np.nan)
outside = np.flatnonzero(~shapely.contains_xy(geometry, target[:, 0], target[:, 1]))
if not len(outside):
return slid
lines = shapely.shortest_line(shapely.points(target[outside]), edge)
nearest = shapely.get_coordinates(lines).reshape(-1, 2, 2)[:, 1]
inward = nearest - target[outside]
length = np.linalg.norm(inward, axis=1)
off_edge = length > 0.0
slid[outside[off_edge]] = (
nearest[off_edge]
+ inward[off_edge] * (distance_tolerance / length[off_edge])[:, None]
)
return slid
def _confine_step(
geometry,
edge,
current,
candidate,
distance_tolerance,
clearance=None,
reach=None,
):
"""Keep moves inside ``geometry``, sliding blocked ones along its edge.
A move whose path leaves ``geometry`` ends at the point of ``edge``
nearest its target, ``distance_tolerance`` inside, when the path there
stays inside and within ``reach``. Other blocked moves stop at their
first crossing, so no move crosses a hole or a gap between parts.
``clearance`` is each turbine's distance from ``edge``, and a move
shorter than it stays inside without a test. ``reach`` holds the initial
positions and the largest distance a turbine may move from them.
"""
tested = np.arange(len(current))
if clearance is not None:
length = np.linalg.norm(candidate - current, axis=1)
# The margin covers rounding in the distance, so a move it lets pass
# is one the containment test would pass.
tested = np.flatnonzero(length >= clearance * (1.0 - 1e-9))
if not len(tested):
return candidate
paths = shapely.linestrings(np.stack((current[tested], candidate[tested]), axis=1))
blocked = tested[~shapely.covers(geometry, paths)]
if not len(blocked):
return candidate
start = current[blocked]
slid = _slid_to_edge(geometry, edge, candidate[blocked], distance_tolerance)
can_slide = ~np.isnan(slid[:, 0])
if reach is not None:
initial, maximum = reach
distance = np.linalg.norm(slid - initial[blocked], axis=1)
can_slide &= distance <= maximum
sliding = np.flatnonzero(can_slide)
if len(sliding):
paths = shapely.linestrings(np.stack((start[sliding], slid[sliding]), axis=1))
sliding = sliding[shapely.covers(geometry, paths)]
candidate[blocked[sliding]] = slid[sliding]
cut = np.ones(len(blocked), dtype=bool)
cut[sliding] = False
blocked, start = blocked[cut], start[cut]
displacement = candidate[blocked] - start
length = np.linalg.norm(displacement, axis=1)
lower = np.zeros(len(blocked))
upper = np.ones(len(blocked))
# Segment containment handles holes and disconnected farms.
active = np.flatnonzero(
(length > 0.0) & ((upper - lower) * length > distance_tolerance)
)
while len(active):
fraction = 0.5 * (lower[active] + upper[active])
point = start[active] + fraction[:, None] * displacement[active]
segments = shapely.linestrings(np.stack((start[active], point), axis=1))
covered = shapely.covers(geometry, segments)
lower[active[covered]] = fraction[covered]
upper[active[~covered]] = fraction[~covered]
remaining = (upper[active] - lower[active]) * length[
active
] > distance_tolerance
active = active[remaining]
moved = length > 0.0
candidate[blocked[moved]] = start[moved] + lower[moved, None] * displacement[moved]
return candidate
def _centering_energy(positions, target, spacing):
"""Return dimensionless centroid-offset energy and a gradient callable."""
if not len(positions):
return 0.0, _zero_gradient(positions)
offset = positions.mean(axis=0) - target
energy = 0.5 * float(offset @ offset) / spacing**2
def gradient():
"""Return the centroid-offset gradient at ``positions``."""
return np.broadcast_to(
offset / (len(positions) * spacing**2), positions.shape
).copy()
return energy, gradient
def _held_off(current, candidate, minimum, distance_tolerance):
"""``candidate`` with moves halved until no pair is closer than ``minimum``.
Each round halves the moves of the turbines in a crowded pair, and a move
halved to ``distance_tolerance`` is dropped. No pair of ``current`` is
closer, so the rounds end. A point along a confined move stays inside
the boundary and the displacement limit, so a halved move does too.
"""
from scipy.spatial import cKDTree # ty: ignore[unresolved-import]
move = candidate - current
length = np.linalg.norm(move, axis=1)
moving = length > 0.0
while np.any(moving):
pairs = cKDTree(candidate).query_pairs(minimum, output_type="ndarray")
differences = candidate[pairs[:, 0]] - candidate[pairs[:, 1]]
crowded = pairs[np.einsum("ij,ij->i", differences, differences) < minimum**2]
crowded = np.unique(crowded)
crowded = crowded[moving[crowded]]
if not len(crowded):
break
move[crowded] *= 0.5
length[crowded] *= 0.5
dropped = crowded[length[crowded] <= distance_tolerance]
move[dropped] = 0.0
moving[dropped] = False
candidate[crowded] = current[crowded] + move[crowded]
return candidate
[docs]
class Relaxation(LayoutModifier):
r"""Relax an existing layout inside its boundary.
``method="repulsion"`` uses the all-pairs potential
.. math:: E_{ij}=\frac{h^2}{d_{ij}^2+\epsilon^2}.
Every pair pushes apart, so the layout expands into the available
geometry. ``method="gaussian"`` uses the shorter-range potential
.. math:: E_{ij}=\exp(-d_{ij}^2/h^2),
plus a penalty for pairs closer than ``length_scale`` and the centroid energy
.. math:: E_c=\frac{\|\bar p-c_B\|^2}{2h^2},
where ``c_B`` is the boundary centroid. Steps follow the per-turbine
normalized negative gradient and accept only energy decreases. A step
that would leave the boundary slides to the nearest point of its edge,
unless that crosses a hole or a gap between parts, where it stops at the
crossing instead. Turbine count and ordering are unchanged. The repulsion
energy visits every turbine pair in cache-sized blocks, so its time and
memory grow quadratically with the turbine count. The Gaussian energy
only visits pairs within four length scales, as farther pairs hold less than
``exp(-16)`` of it, and gradients are only computed for accepted steps.
Parameters
----------
method:
``"repulsion"`` for inverse-square spreading, ``"gaussian"`` for
local pair repulsion.
length_scale:
The length ``h`` of the potentials and the compact-penalty threshold,
in meters. Defaults to the square root of boundary area per turbine.
max_iterations:
Maximum number of energy-descent steps.
step:
Maximum movement per iteration as a fraction of ``length_scale``.
tolerance:
Stop once the energy falls by at most this fraction over five
accepted steps.
max_displacement:
Optional maximum total movement from each initial position, in meters.
``None`` leaves displacement unconstrained.
minimum_spacing:
Optional distance in meters no step brings a pair closer than. The
moves of a pair that would come closer are halved until it does not.
A layout filled at a minimum spacing holds its relaxation to it.
seed:
Seed used only to separate exactly coincident positions.
"""
def __init__(
self,
*,
method: typing.Literal["repulsion", "gaussian"] = "repulsion",
length_scale: float | None = None,
max_iterations: int = 100,
step: float = 0.1,
tolerance: float = 1e-3,
max_displacement: float | None = None,
minimum_spacing: float | None = None,
seed: int | np.random.Generator | None = None,
) -> None:
if method not in {"repulsion", "gaussian"}:
raise ValueError("method must be 'repulsion' or 'gaussian'")
self.method = method
self.length_scale = (
None
if length_scale is None
else _validate_positive("length_scale", length_scale)
)
self.max_iterations = _validate_integer("max_iterations", max_iterations)
self.step = _validate_positive("step", step)
if self.step > 1:
raise ValueError("step must be in (0, 1]")
self.tolerance = _validate_positive("tolerance", tolerance)
self.max_displacement = (
None
if max_displacement is None
else _validate_nonnegative("max_displacement", max_displacement)
)
self.minimum_spacing = (
None
if minimum_spacing is None
else _validate_positive("minimum_spacing", minimum_spacing)
)
self.seed = pin_seed(seed)
def _held_apart(self, spacing: float) -> "Relaxation":
held = copy(self)
held.minimum_spacing = max(self.minimum_spacing or 0.0, spacing)
return held
def _spacing(self, layout: Farm) -> float:
if self.length_scale is not None:
return self.length_scale
boundary = layout.boundary
if boundary is None:
raise ValueError("an automatic length_scale requires a farm boundary")
return float(np.sqrt(boundary.area / max(layout.n_turbines, 1)))
[docs]
def apply(self, layout: Farm, rng: np.random.Generator, seed: int) -> Farm:
"""Return a new farm with positions relaxed inside the boundary.
``rng`` only separates exactly coincident positions. Records the
energies, largest displacement and convergence in
``metadata["modifiers"]``. The string guide is dropped when any turbine
moves. Raises ``ValueError`` without a boundary, when a position
starts outside it, or when a pair starts closer than
``minimum_spacing``.
"""
if layout.boundary is None:
raise ValueError("relaxation requires a farm boundary")
if not np.all(layout.boundary._contains(layout.positions)):
raise ValueError(
"relaxation requires all initial positions inside the boundary"
)
# The margin covers rounding in the distances a fill measured.
if self.minimum_spacing is not None and layout.minimum_spacing() < (
self.minimum_spacing * (1.0 - 1e-9)
):
raise ValueError("layout does not satisfy the requested minimum_spacing")
initial = np.asarray(layout.positions, dtype=float)
positions = initial.copy()
spacing = self._spacing(layout)
distance_tolerance = _DISTANCE_TOLERANCE * spacing
centering_target = np.asarray(layout.boundary.centroid)
if self.method == "repulsion":
softening = distance_tolerance
def evaluate(points: np.ndarray) -> tuple[float, Callable[[], np.ndarray]]:
"""Return the inverse-square energy and its uphill gradient callable."""
return _repulsion_energy(points, spacing, rng, softening)
else:
def evaluate(points: np.ndarray) -> tuple[float, Callable[[], np.ndarray]]:
"""Return the Gaussian, spacing and centroid energy and gradient callable."""
pair_energy, pair_gradient = _gaussian_energy(points, spacing, rng)
center_energy, center_gradient = _centering_energy(
points, centering_target, spacing
)
def total_gradient() -> np.ndarray:
"""Return the summed uphill gradient at ``points``."""
return pair_gradient() + center_gradient()
return pair_energy + center_energy, total_gradient
# Gradients are only computed for accepted steps.
initial_energy, gradient_at = evaluate(positions)
gradient = gradient_at()
energy = initial_energy
converged = layout.n_turbines < 2 or self.max_displacement == 0.0
iteration = 0
# Every trial step tests its moves against the boundary.
geometry = layout.boundary.geometry
edge = geometry.boundary
shapely.prepare(geometry)
shapely.prepare(edge)
reach = (
None if self.max_displacement is None else (initial, self.max_displacement)
)
# The line search starts at the step last accepted, and at twice it
# once a step passes at its first try, never above the full step. A
# step that had to halve mostly halves again, so starting above it
# wastes an evaluation.
start_fraction = 1.0
history = [energy]
for iteration in range(1, self.max_iterations + 1):
if converged:
iteration = 0
break
# The gradient at ``positions`` came with its accepted energy.
gradient_norm = np.linalg.norm(gradient, axis=1)
if float(gradient_norm.max(initial=0.0)) == 0.0:
converged = True
break
direction = np.divide(
gradient,
gradient_norm[:, None],
out=np.zeros_like(gradient),
where=gradient_norm[:, None] > 0.0,
)
displacement = -self.step * spacing * direction
clearance = shapely.distance(edge, shapely.points(positions))
accepted = False
fraction = start_fraction
movement = 0.0
while (
fraction * float(np.linalg.norm(displacement, axis=1).max(initial=0.0))
> distance_tolerance
):
candidate = positions + fraction * displacement
candidate = _limit_displacement(
candidate, initial, self.max_displacement
)
candidate = _confine_step(
geometry,
edge,
positions,
candidate,
distance_tolerance,
clearance,
reach,
)
if self.minimum_spacing is not None:
candidate = _held_off(
positions, candidate, self.minimum_spacing, distance_tolerance
)
candidate_energy, candidate_gradient = evaluate(candidate)
if candidate_energy <= energy:
movement = float(
np.linalg.norm(candidate - positions, axis=1).max(initial=0.0)
)
positions = candidate
energy, gradient_at = candidate_energy, candidate_gradient
gradient = gradient_at()
accepted = True
break
fraction *= 0.5
if fraction == start_fraction:
start_fraction = min(1.0, 2.0 * fraction)
else:
start_fraction = fraction
if not accepted or movement <= distance_tolerance:
converged = True
break
history.append(energy)
if len(history) > _STALL_STEPS:
earlier = history[-_STALL_STEPS - 1]
if earlier - energy <= self.tolerance * abs(earlier):
converged = True
break
displacements = np.linalg.norm(positions - initial, axis=1)
metadata = dict(layout.metadata)
modifiers = list(metadata.get("modifiers", []))
record = {
"type": "relaxation",
"method": self.method,
"seed": int(seed),
"length_scale": spacing,
"iterations": iteration,
"maximum_iterations": self.max_iterations,
"step": self.step,
"tolerance": self.tolerance,
"max_displacement": self.max_displacement,
"minimum_spacing": self.minimum_spacing,
"maximum_displacement": float(displacements.max(initial=0.0)),
"initial_energy": initial_energy,
"final_energy": energy,
"converged": converged,
}
if self.method == "gaussian":
record["centering_target"] = centering_target.tolist()
modifiers.append(record)
metadata["modifiers"] = modifiers
moved = bool(np.any(displacements > 0.0))
return Farm(
positions,
layout.boundary,
metadata=metadata,
guide=None if moved else layout.guide,
resource=layout.resource,
turbines=layout.turbines,
surface=layout.surface,
)
