Source code for oowl.layouts.uniform

"""Uniformly random turbine layouts."""

from __future__ import annotations

import numpy as np
import shapely

from .base import Farm, LayoutAlgorithm

__all__ = ["Uniform"]


def _triangulate(geometry) -> np.ndarray:
    """Return the positive-area triangles of a constrained triangulation.

    Area-weighted triangle choice then a uniform point inside is exactly
    uniform over ``geometry``, holes included.
    """
    triangles = shapely.get_parts(shapely.constrained_delaunay_triangles(geometry))
    return triangles[shapely.area(triangles) > 0]


def _uniform_points(geometry, count: int, rng: np.random.Generator) -> np.ndarray:
    """``count`` independent uniform points strictly inside ``geometry``.

    Raises ``ValueError`` for a geometry without finite positive area.
    """
    triangles = _triangulate(geometry)
    if not len(triangles):
        raise ValueError("Uniform requires a boundary with finite positive area")

    areas = shapely.area(triangles)
    total_area = float(areas.sum())
    if not np.isfinite(total_area) or total_area <= 0:
        raise ValueError("Uniform requires a boundary with finite positive area")

    selected_triangles = rng.choice(len(triangles), size=count, p=areas / total_area)
    # Every triangle's closed ring holds four points, its corners and the
    # first again. Gathered one drawn point at a time, they took a third of
    # the time a mixture spent drawing its cells.
    corners = shapely.get_coordinates(shapely.get_exterior_ring(triangles))
    vertices = corners.reshape(len(triangles), 4, 2)[selected_triangles, :3]
    # Uniform barycentric sampling of a triangle, from Osada et al. (2002),
    # https://doi.org/10.1145/571647.571648. Clamping the draws off 0 and 1
    # keeps points strictly inside, off the triangle's edges.
    draws = rng.random((count, 2))
    open_zero = np.nextafter(0.0, 1.0)
    first = np.sqrt(np.maximum(draws[:, 0], open_zero))
    second = np.clip(draws[:, 1], open_zero, np.nextafter(1.0, 0.0))
    return (
        (1.0 - first)[:, None] * vertices[:, 0]
        + (first * (1.0 - second))[:, None] * vertices[:, 1]
        + (first * second)[:, None] * vertices[:, 2]
    )




class Uniform(LayoutAlgorithm):
    """Scatter turbines independently and uniformly across the boundary.

    Every position is an independent draw with no lattice and no minimum
    separation, so pairs may fall arbitrarily close together. A
    ``minimum_spacing`` therefore sizes the farm by its closest pair, which
    one unlucky draw sets, so the area or count it gives varies widely
    between seeds. This is complete spatial randomness, the null model behind
    :func:`~oowl.metrics.clark_evans`, :func:`~oowl.metrics.ripley_k` and
    :func:`~oowl.metrics.quadrat_p_value`. Unrelated to the wind resource
    :class:`oowl.resources.UniformDirection`.

    Parameters
    ----------
    n_turbines:
        Exact number of turbines to generate.
    density:
        The one density setting the count, ``{"minimum_spacing": meters}``
        for the most turbines with no pair closer, ``{"mean_spacing":
        meters}`` for the most whose :meth:`~oowl.Farm.mean_spacing` is at
        least it, ``{"maximum_spacing": meters}`` for the fewest leaving none
        farther from a neighbor, or ``{"area_per_turbine": square meters}``.
    seed:
        Random seed controlling the draws.
    """

    def __init__(
        self,
        n_turbines: int | None = None,
        *,
        density: dict[str, float] | None = None,
        seed: int | np.random.Generator | None = None,
    ) -> None:
        super().__init__(
            n_turbines,
            density=density,
            seed=seed,
        )

    def _fill_estimate(self, boundary, spacing: float, mean: bool = False) -> int:
        """The count whose mean nearest-neighbor distance is about ``spacing``.

        Independent uniform points of density ``n / area`` lie a mean of
        ``0.5 * sqrt(area / n)`` from their nearest neighbor. The edge leaves
        points near it fewer neighbors, so a farm holds a few more.
        """
        return max(2, int(float(boundary.area) / (4.0 * spacing**2)))

    def _generate_single(
        self,
        boundary,
        rng: np.random.Generator,
        seed: int,
    ) -> Farm:
        positions = _uniform_points(boundary.geometry, self._required_n, rng)
        return self._single_result(
            positions,
            boundary,
            seed,
            {"method": "uniform"},
        )