"""Aligned parallelogram-grid layout."""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
from .._validation import _validate_finite, _validate_positive
from .candidate import _CandidateLayout, _ray_counter
if TYPE_CHECKING: # pragma: no cover
from ..boundaries.base import Boundary
from .fitters import SelectedCandidateFitter
__all__ = ["Grid"]
[docs]
class Grid(_CandidateLayout):
"""Occupy selected intersections of a regular parallelogram grid.
``spacing_ratio`` is the length of the second grid vector divided by the
first. ``lattice_angle`` is the angle in degrees between them, and
``rotation`` orients the first vector globally. The fitter finds the
largest aligned grid that offers at least ``n_turbines`` valid
intersections, then selects an exact, well-distributed subset.
"""
def __init__(
self,
n_turbines: int | None = None,
*,
density: dict[str, float] | None = None,
spacing_ratio: float = 1.0,
lattice_angle: float = 90.0,
rotation: float = 0.0,
fitter: SelectedCandidateFitter | None = None,
seed: int | np.random.Generator | None = None,
):
self.spacing_ratio = _validate_positive("spacing_ratio", spacing_ratio)
self.lattice_angle = _validate_finite("lattice_angle", lattice_angle)
if not 0 < self.lattice_angle < 180:
raise ValueError("lattice_angle must be between 0 and 180 degrees")
self.rotation = _validate_finite("rotation", rotation)
super().__init__(
n_turbines,
density=density,
fitter=fitter,
seed=seed,
)
def _basis(self, scale: float) -> np.ndarray:
angle = np.deg2rad(self.lattice_angle)
basis = scale * np.array(
[
[1.0, 0.0],
[
self.spacing_ratio * np.cos(angle),
self.spacing_ratio * np.sin(angle),
],
]
)
if self.rotation:
theta = np.deg2rad(self.rotation)
cosine, sine = np.cos(theta), np.sin(theta)
rotation = np.array([[cosine, -sine], [sine, cosine]])
basis = basis @ rotation.T
return basis
[docs]
def estimate_spacing(self, boundary: Boundary) -> float:
"""Return a starting length of the first grid vector in meters.
Assumes one parallelogram cell per turbine. The estimate only seeds the
fitter's search and is not the final spacing.
"""
cell_factor = self.spacing_ratio * np.sin(np.deg2rad(self.lattice_angle))
return float(np.sqrt(boundary.area / (self._required_n * cell_factor)))
def _detailed_fit(self, boundary, positions, fit):
# Markov turbine removal finds grid rows from the basis.
return positions, {**fit, "grid_basis": self._basis(fit["scale"]).tolist()}
def _lattice(self, boundary, scale, phases):
"""The basis, origin, phase offsets and index ranges at ``scale``.
Each phase's inclusive index ranges cover the boundary's box.
"""
basis = self._basis(scale)
inverse = np.linalg.inv(basis)
center = np.asarray(boundary.centroid, dtype=float)
minx, miny, maxx, maxy = boundary.bounds
# Lattice offsets of the bounding-box corners from the grid origin.
corners = (
np.array(
[
[minx, miny],
[minx, maxy],
[maxx, miny],
[maxx, maxy],
]
)
- center
)
offsets = np.asarray(phases, dtype=float).reshape(-1, 2)
# Elementwise rather than matrix products, because BLAS kernels vary
# with the row count and a product would make a lattice point's last
# bits depend on how many other points were generated beside it.
translations = offsets[:, :1] * basis[0] + offsets[:, 1:] * basis[1]
coefficients = (corners[None, :, :] - translations[:, None, :]) @ inverse
lower = np.floor(coefficients.min(axis=1)).astype(int) - 1
upper = np.ceil(coefficients.max(axis=1)).astype(int) + 1
return basis, center, offsets, translations, lower, upper
def _candidate_bounds(self, boundary, scale, phases, max_candidates):
if scale <= 0:
return None
_, center, offsets, _, lower, upper = self._lattice(boundary, scale, phases)
shape = upper - lower + 1
# A larger scale narrows each index range, give or take one index
# each side for rounding, so no field above ``scale`` passes the cap.
if np.any((shape[:, 0] + 2) * (shape[:, 1] + 2) > max_candidates):
return None
unit = self._basis(1.0)
# Each phase's index range, its first index varying fastest.
cells = shape[:, 0] * shape[:, 1]
phase = np.repeat(np.arange(len(offsets)), cells)
cell = np.arange(len(phase)) - np.repeat(np.cumsum(cells) - cells, cells)
along = (cell % shape[phase, 0] + lower[phase, 0]) + offsets[phase, 0]
across = (cell // shape[phase, 0] + lower[phase, 1]) + offsets[phase, 1]
units = np.multiply.outer(along, unit[0]) + np.multiply.outer(across, unit[1])
return _ray_counter(boundary, center, units, phase, len(offsets), scale)
def _candidate_fields(
self, boundary, scale, phases, max_candidates=None, min_candidates=None
):
if scale <= 0:
return [np.empty((0, 2)) for _ in phases]
basis, center, offsets, translations, lower, upper = self._lattice(
boundary, scale, phases
)
phase, second, first_low, counts = _clipped_rows(
boundary, basis, center, translations, lower, upper
)
# Candidates are only built inside the boundary's rotated rectangle,
# so a thin or diagonal boundary costs what it holds, not its box. Rows
# still long and sparse there, as across a thin L, are also cut at the
# boundary's edges where that costs less than testing their points.
runs = _edge_clipped_runs(
boundary, basis, center, translations, phase, second, first_low, counts
)
if runs is not None:
phase, second, first_low, counts = runs
cells = np.bincount(phase, weights=counts, minlength=len(offsets)).astype(int)
within = np.ones(len(offsets), dtype=bool)
if max_candidates is not None:
within = cells <= max_candidates
if not np.any(within):
return [None for _ in phases]
fields = [np.empty((0, 2)) if is_within else None for is_within in within]
# A phase's rows bound its candidates, so a phase that cannot reach
# ``min_candidates`` skips the containment query.
computed = (
within if min_candidates is None else within & (cells >= min_candidates)
)
kept = computed[phase] & (counts > 0)
if not np.any(kept):
return fields
lengths = counts[kept]
starts = np.cumsum(lengths) - lengths
first = np.repeat(first_low[kept], lengths) + (
np.arange(int(lengths.sum())) - np.repeat(starts, lengths)
)
point_phase = np.repeat(phase[kept], lengths)
# The same elementwise sums as a whole lattice, so each point keeps its
# last bits wherever it was built.
lattice = np.multiply.outer(first, basis[0]) + np.multiply.outer(
np.repeat(second[kept], lengths), basis[1]
)
points = lattice + translations[point_phase]
points += center
inside = boundary._contains(points)
# Rows run phase by phase, each in index order, as a whole lattice does,
# so each computed phase holds one run of the points.
stops = np.cumsum(cells[computed])
for index, start, stop in zip(
np.flatnonzero(computed), stops - cells[computed], stops
):
fields[index] = points[start:stop][inside[start:stop]]
return fields
def _clipped_rows(boundary, basis, center, translations, lower, upper):
"""Each phase's lattice rows clipped to the boundary's rotated rectangle.
Returns per row its phase, second index, first first index and count of
points, phase by phase and each in index order. A point left out lies
outside the rectangle and so outside the boundary.
"""
spans = upper[:, 1] - lower[:, 1] + 1
phase = np.repeat(np.arange(len(translations)), spans)
second = (
np.arange(len(phase)) - np.repeat(np.cumsum(spans) - spans, spans)
) + lower[phase, 1]
origins = center + translations[phase] + second[:, None] * basis[1]
normals, limits = boundary._enclosing_halfplanes
along = normals @ basis[0]
room = limits[None, :] - origins @ normals.T
low = lower[phase, 0].astype(float)
high = upper[phase, 0].astype(float)
rising, falling = along > 0, along < 0
if rising.any():
high = np.minimum(high, np.floor((room[:, rising] / along[rising]).min(axis=1)))
if falling.any():
low = np.maximum(low, np.ceil((room[:, falling] / along[falling]).max(axis=1)))
# A row parallel to a side lies wholly inside or outside it.
level = ~(rising | falling)
if level.any():
high[(room[:, level] < 0).any(axis=1)] = -np.inf
counts = np.maximum(high - low + 1, 0).astype(int)
return phase, second, low.astype(int), counts
#: Cost of cutting one row at one edge, as a share of testing one point,
#: which decides whether rows are cut at the boundary's edges.
_EDGE_PAIR_COST = 0.25
def _edge_clipped_runs(
boundary, basis, center, translations, phase, second, low, counts
):
"""The runs of each row that may lie in the boundary, cut at its edges.
Each edge reaching a row, widened by a margin, marks a stretch whose
points may lie inside. Between two stretches the row meets no edge, so it
lies wholly inside or outside, which one containment test at its middle
tells, and an inside gap joins its stretches. Returns runs as
:func:`_clipped_rows` returns rows, several to a row, each a superset of
the row's points inside, so containment still decides every point.
``None`` when cutting would cost more than testing the rows' points.
"""
rows = np.flatnonzero(counts)
empty = np.empty(0, dtype=int)
if not len(rows):
return empty, empty, empty, empty
starts, ends = boundary._edges
along = basis[0]
length = float(np.hypot(*along))
across = np.array([-along[1], along[0]]) / length
origins = center + translations[phase[rows]] + second[rows, None] * basis[1]
minx, miny, maxx, maxy = boundary.bounds
margin = 1e-9 * float(np.hypot(maxx - minx, maxy - miny)) + 1e-6
# The rows are parallel, so an edge reaches those whose line lies within
# its extent across them, a run of the rows sorted by that height.
height = origins @ across
by_height = np.argsort(height, kind="stable")
sorted_height = height[by_height]
start_height, end_height = starts @ across, ends @ across
first_row = np.searchsorted(
sorted_height, np.minimum(start_height, end_height) - margin, side="left"
)
last_row = np.searchsorted(
sorted_height, np.maximum(start_height, end_height) + margin, side="right"
)
reached = last_row - first_row
crossings = int(reached.sum())
if crossings * _EDGE_PAIR_COST >= counts.sum():
return None
if not crossings:
return empty, empty, empty, empty
edge = np.repeat(np.arange(len(starts)), reached)
row = by_height[
np.arange(crossings)
- np.repeat(np.cumsum(reached) - reached, reached)
+ np.repeat(first_row, reached)
]
rise_from = start_height[edge] - height[row]
rise_to = end_height[edge] - height[row]
rise = rise_to - rise_from
# The share of each edge within the margin of the line, all of one along it.
level = rise == 0.0
safe = np.where(level, 1.0, rise)
first = np.where(level, 0.0, (-margin - rise_from) / safe)
last = np.where(level, 1.0, (margin - rise_from) / safe)
enter = np.clip(np.minimum(first, last), 0.0, 1.0)
leave = np.clip(np.maximum(first, last), 0.0, 1.0)
# Positions along the row in lattice steps from its origin.
offset = (origins @ along)[row] / length**2
begin = (starts[edge] @ along) / length**2 - offset
finish = (ends[edge] @ along) / length**2 - offset
slack = 2.0 * margin / length
one = begin + enter * (finish - begin)
two = begin + leave * (finish - begin)
lows = np.minimum(one, two) - slack
highs = np.maximum(one, two) + slack
# An edge whose ends both lie past the margin crosses the line cleanly, at
# one point, so it turns the row from outside to inside or back. One with
# an end within the margin, as at a vertex, leaves the row's side unknown.
clean = ((rise_from < -margin) & (rise_to > margin)) | (
(rise_from > margin) & (rise_to < -margin)
)
order = np.lexsort((lows, row))
row, lows, highs, clean = row[order], lows[order], highs[order], clean[order]
# A stretch starting past every one before it in its row opens a block. The
# rows are kept apart by an offset, which only decides where blocks open.
apart = float(highs.max() - lows.min()) + 1.0
reached = np.maximum.accumulate(highs + row * apart) - row * apart
opens = np.ones(len(row), dtype=bool)
opens[1:] = (row[1:] != row[:-1]) | (lows[1:] > reached[:-1])
blocks = np.flatnonzero(opens)
block_row, block_low = row[blocks], lows[blocks]
block_high = np.maximum.reduceat(highs, blocks)
flips = np.add.reduceat(clean.astype(int), blocks)
doubtful = np.maximum.reduceat((~clean).astype(int), blocks).astype(bool)
# A gap between two blocks of a row meets no edge, so it lies wholly
# inside or outside. The row starts outside, and each block turns it by
# its clean crossings. After a doubtful block the gap's middle is tested,
# and the turns count on from there.
index = np.arange(len(blocks))
first_block = np.maximum.accumulate(
np.where(np.append(True, block_row[1:] != block_row[:-1]), index, 0)
)
turned = np.cumsum(flips)
turned = turned - (turned[first_block] - flips[first_block])
anchor = np.maximum.accumulate(np.where(doubtful, index, -1))
anchored = anchor >= first_block
same = np.append(block_row[1:] == block_row[:-1], False)
inside = np.zeros(len(blocks), dtype=bool)
tested = same & anchored & (anchor == index)
if tested.any():
middles = 0.5 * (block_high[tested] + block_low[np.flatnonzero(tested) + 1])
points = origins[block_row[tested]] + middles[:, None] * along
inside[tested] = boundary._contains(points)
counted = same & ~tested
base = np.where(anchored, inside[np.maximum(anchor, 0)], False)
since = np.where(anchored, turned - turned[np.maximum(anchor, 0)], turned)
inside[counted] = base[counted] ^ (since[counted] % 2 == 1)
joined = inside[:-1]
runs = np.flatnonzero(np.append(True, ~joined))
run_row = block_row[runs]
run_low = block_low[runs]
run_high = np.maximum.reduceat(block_high, runs)
source = rows[run_row]
start = np.maximum(np.ceil(run_low), low[source]).astype(int)
stop = np.minimum(np.floor(run_high), low[source] + counts[source] - 1).astype(int)
return phase[source], second[source], start, np.maximum(stop - start + 1, 0)
