Source code for structure_factor.spatial_windows

"""Collection of classes representing observation windows (box, ball, etc).

- :py:class:`~structure_factor.spatial_windows.BallWindow`: Ball window object.
- :py:class:`~structure_factor.spatial_windows.BoxWindow`: Box window object.

.. note::

    **Typical usage**

    - :py:class:`~structure_factor.point_pattern.PointPattern` has a :py:attr:`~structure_factor.point_pattern.PointPattern.window` argument/attribute.
"""

from abc import ABCMeta, abstractmethod

import numpy as np
import scipy as sp
from matplotlib import pyplot as plt
from matplotlib.patches import Circle, Rectangle
from rpy2 import robjects
from spatstat_interface.interface import SpatstatInterface

from structure_factor.utils import get_random_number_generator


[docs]class AbstractSpatialWindow(metaclass=ABCMeta): r"""Encapsulate the notion of spatial window in :math:`\mathbb{R}^d`.""" @property @abstractmethod def dimension(self): """Return the ambient dimension of the corresponding window.""" @property @abstractmethod def volume(self): """Compute the volume of the corresponding window.""" @abstractmethod def __contains__(self, point): r"""Return True if :math:`d` dimensional vector ``point`` falls inside the corresponding :math:`d` dimensional window, otherwise return False. Args: point (numpy.ndarray): :math:`d` dimensional vector to test. """
[docs] def indicator_function(self, points): r"""Return the indicator function of the corresponding window evaluated at each of the :math:`n` ``points``. Args: points (numpy.ndarray): Vector of size :math:`d` or array of size :math:`n \times d` containing the point(s) to be tested. Returns: bool or numpy.ndarray: - If :math:`n=1`, bool. - If :math:`n>1`, :math:`n` dimensional boolean array. """ if points.ndim == 1 and points.size == self.dimension: return points in self return np.apply_along_axis(self.__contains__, axis=1, arr=points)
[docs] @abstractmethod def rand(self, n=1, seed=None): r"""Generate `n` points uniformly at random in the corresponding spatial window. Args: n (int, optional): Number of points. Defaults to 1. seed (int or np.random.Generator, optional): Defaults to None. Returns: numpy.ndarray: - If :math:`n=1`, :math:`d` dimensional vector. - If :math:`n>1`, :math:`n \times d` array containing the points. """
[docs]class BallWindow(AbstractSpatialWindow): r"""Create a :math:`d` dimensional ball window :math:`B(c, r)`, where :math:`c \in \mathbb{R}^d` and :math:`r>0`. .. todo:: list attributes Example: .. plot:: code/spatial_window/ball_window.py :include-source: True :align: center .. seealso:: - :py:mod:`~structure_factor.point_pattern` - :py:class:`~structure_factor.spatial_windows.BoxWindow` """
[docs] def __init__(self, center, radius=1.0): """Initialize a :math:`d` dimensional ball window :math:`B(c, r)` from the prescribed ``center`` and ``radius``. Args: center (numpy.ndarray): Center :math:`c` of the ball. radius (float, optional): Radius :math:`r > 0` of the ball. Defaults to 1.0. """ center = np.asarray(center) if not center.ndim == 1: raise ValueError("center must be 1D numpy.ndarray") if not radius > 0: raise ValueError("radius must be positive") self.center = center self.radius = float(radius)
@property def dimension(self): return len(self.center) @property def surface(self): d, r = self.dimension, self.radius if d == 1: return 0.0 if d == 2: return 2 * np.pi * r if d == 3: return 4 * np.pi * r ** 2 return 2 * np.pi ** (d / 2) * r ** (d - 1) / sp.special.gamma(d / 2) @property def volume(self): d, r = self.dimension, self.radius if d == 1: return 2 * r if d == 2: return np.pi * r ** 2 if d == 3: return 4 / 3 * np.pi * r ** 3 return np.pi ** (d / 2) * r ** d / sp.special.gamma(d / 2 + 1) def __contains__(self, point): point = np.asarray(point) assert point.ndim == 1 and point.size == self.dimension return self.indicator_function(point)
[docs] def indicator_function(self, points): return np.linalg.norm(points - self.center, axis=-1) <= self.radius
[docs] def rand(self, n=1, seed=None): # Method of dropped coordinates # Efficiently sampling vectors and coordinates from the n-sphere and n-ball # Voelker, Aaron and Gosmann, Jan and Stewart, Terrence # doi: 10.13140/RG.2.2.15829.01767/1 rng = get_random_number_generator(seed) d = self.dimension points = rng.standard_normal(size=(n, d + 2)) points /= np.linalg.norm(points, axis=-1, keepdims=True) idx = 0 if n == 1 else slice(0, n) return self.center + self.radius * points[idx, :d]
[docs] def to_spatstat_owin(self, **params): """Convert the object to a ``spatstat.geom.disc`` R object of type ``disc``, which is a subtype of ``owin``. Args: params (dict): Optional keyword arguments passed to ``spatstat.geom.disc``. Returns: spatstat.geom.disc: R object. .. seealso:: - `https://rdocumentation.org/packages/spatstat.geom/versions/2.2-0/topics/disc <https://rdocumentation.org/packages/spatstat.geom/versions/2.2-0/topics/disc>`_ """ spatstat = SpatstatInterface(update=False) spatstat.import_package("geom", update=False) r = self.radius c = robjects.vectors.FloatVector(self.center) return spatstat.geom.disc(radius=r, centre=c, **params)
[docs] def plot(self, axis=None, **kwargs): """Display the window on matplotlib `axis`. Args: axis (plt.Axes, optional): Support axis of the plot. Defaults to None. Keyword Args: kwargs (dict): Keyword arguments of ``matplotlib.patches.Circle`` with default ``fill=False``. Returns: plt.Axes: Plot axis. """ if self.dimension != 2: raise NotImplementedError("Method implemented only for 2D window") if axis is None: fig, axis = plt.subplots(figsize=(5, 5)) kwargs.setdefault("fill", False) circle = Circle(self.center, self.radius, **kwargs) axis.add_patch(circle) return axis
[docs]class UnitBallWindow(BallWindow): r"""Create a d-dimensional unit ball window :math:`B(c, r=1)`, where :math:`c \in \mathbb{R}^d`. .. note:: ``UnitBallWindow(center) = BallWindow(center, radius=1.0)`` """
[docs] def __init__(self, center): """Initialize a :math:`d` dimensional unit ball window :math:`B(c, r=1)` from the prescribed ``center``. Args: center (numpy.ndarray, optional): Center :math:`c` of the ball. """ super().__init__(center, radius=1.0)
[docs]class BoxWindow(AbstractSpatialWindow): r"""Create a :math:`d` dimensional box window :math:`\prod_{i=1}^{d} [a_i, b_i]`. .. todo:: list attributes Example: .. plot:: code/spatial_window/box_window.py :include-source: True :align: center .. seealso:: - :py:mod:`~structure_factor.point_pattern` - :py:class:`~structure_factor.spatial_windows.BoxWindow` """
[docs] def __init__(self, bounds): r"""Initialize :math:`d` dimensional unit box window the prescibed ``bounds[i, :]`` :math:`=[a_i, b_i]`. Args: bounds (numpy.ndarray): :math:`d \times 2` array describing the bounds of the box. """ bounds = np.atleast_2d(bounds) if bounds.ndim != 2 or bounds.shape[1] != 2: raise ValueError("bounds must be d x 2 numpy.ndarray") if np.any(np.diff(bounds, axis=-1) <= 0): raise ValueError("all bounds [a_i, b_i] must satisfy a_i < b_i") # use transpose to facilitate operations (unpacking, diff, rand, etc) self._bounds = np.transpose(bounds)
@property def bounds(self): r"""Return the bounds decribing the BoxWindow. ``bounds[i, :]`` :math:`=[a_i, b_i]`. """ return np.transpose(self._bounds) @property def dimension(self): return self._bounds.shape[1] @property def volume(self): return np.prod(np.diff(self._bounds, axis=0)) def __contains__(self, point): point = np.asarray(point) assert point.ndim == 1 and point.size == self.dimension return self.indicator_function(point)
[docs] def indicator_function(self, points): a, b = self._bounds return np.logical_and( np.all(a <= points, axis=-1), np.all(points <= b, axis=-1) )
[docs] def rand(self, n=1, seed=None): rng = get_random_number_generator(seed) a, b = self._bounds d = self.dimension return rng.uniform(a, b, size=(d,) if n == 1 else (n, d))
[docs] def to_spatstat_owin(self, **params): """Convert the object to a ``spatstat.geom.owin`` R object of type ``owin``. Args: params (dict): Optional keyword arguments passed to ``spatstat.geom.owin``. Returns: spatstat.geom.owin: R object. .. seealso:: - `https://rdocumentation.org/packages/spatstat.geom/versions/2.2-0/topics/owin <https://rdocumentation.org/packages/spatstat.geom/versions/2.2-0/topics/owin>`_ """ if self.dimension != 2: raise NotImplementedError("spatstat only handles 2D windows") spatstat = SpatstatInterface(update=False) spatstat.import_package("geom", update=False) a, b = self.bounds x = robjects.vectors.FloatVector(a) y = robjects.vectors.FloatVector(b) return spatstat.geom.owin(xrange=x, yrange=y, **params)
[docs] def plot(self, axis=None, **kwargs): """Display the window on matplotlib `axis`. Args: axis (plt.Axes, optional): Support axis of the plot. Defaults to None. Keyword Args: kwargs (dict): Keyword arguments of ``matplotlib.patches.Rectangle`` with default ``fill=False``. Returns: plt.Axes: Plot axis. """ if self.dimension != 2: raise NotImplementedError("Method implemented only for 2D window") if axis is None: fig, axis = plt.subplots(figsize=(5, 5)) kwargs.setdefault("fill", False) xy = self._bounds[0] width, height = np.diff(self._bounds, axis=0).ravel() rectangle = Rectangle(xy, width, height, **kwargs) axis.add_patch(rectangle) return axis
[docs]class UnitBoxWindow(BoxWindow): r"""Create a :math:`d` dimensional unit box window :math:`\prod_{i=1}^{d} [c_i - \frac{1}{2}, c_i + \frac{1}{2}]` where :math:`c \in \mathbb{R}^d`."""
[docs] def __init__(self, center): r"""Initialize a :math:`d` dimensional unit box window :math:`\prod_{i=1}^{d} [c_i - \frac{1}{2}, c_i + \frac{1}{2}]`, i.e., a box window with length equal to 1 and prescribed ``center``, such that :math:`c_i=` ``center[i]``. Args: center (numpy.ndarray): Center :math:`c` of the box. """ if np.ndim(center) != 1: raise ValueError("center must be 1D array.") bounds = np.add.outer(center, [-0.5, 0.5]) super().__init__(bounds)
[docs]def check_cubic_window(window): """Check whether ``window`` is cubic. Args: window (:py:class:`~structure_factor.spatial_windows.BoxWindow`): Window to be checked. Example: .. literalinclude:: code/spatial_window/check_cubic_window.py :language: python """ if not isinstance(window, BoxWindow): raise TypeError("window must be an instance of BoxWindow.") lengths = np.diff(window.bounds, axis=1) if np.any(lengths != lengths[0]): raise ValueError("window should be a 'cubic' BoxWindow.") return None
[docs]def check_centered_window(window): """Check whether ``window`` is centered at the origin. Args: window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): Window to be checked. Example: .. literalinclude:: code/spatial_window/check_centered_window.py :language: python """ if isinstance(window, BoxWindow): if np.sum(window.bounds) != 0: raise ValueError( "window is not centered at the origin. Hint: use a centered window." ) if isinstance(window, BallWindow): if any(window.center != 0): raise ValueError( "window is not centered at the origin. Hint: use a centered window." ) return None
[docs]def subwindow_parameter_max(window, subwindow_type="BoxWindow"): """ Return the parameter i.e., lengthside (resp. radius) of the largest cubic (resp. ball) subwindow of ``window``. Args: window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): BoxWindow or BallWindow centered at the origin. subwindow_type (str, optional): Type of the subwindow ("BoxWindow" or "BallWindow"). Defaults to "BoxWindow". Returns: float : Parameter of the largest subwindow of ``window``. Example: .. plot:: code/spatial_window/subwindow_param_max.py :include-source: True :align: center """ if subwindow_type not in ["BoxWindow", "BallWindow"]: raise ValueError( "The available subwindow types are BallWindow or BoxWindow. Hint: the parameter corresponding to the window type must be 'BallWindow' or 'BoxWindow'. " ) check_centered_window(window) # window parameter if isinstance(window, BallWindow): if subwindow_type == "BallWindow": param_max = window.radius else: param_max = window.radius * 2 / np.sqrt(2) # length side of the BoxWindow elif isinstance(window, BoxWindow): if subwindow_type == "BallWindow": param_max = np.min(np.diff(window.bounds)) / 2 else: param_max = np.min(np.diff(window.bounds)) return param_max