Source code for psimpy.sampler.latin

import numpy as np
from scipy.spatial.distance import pdist
from typing import Union
from beartype import beartype
from psimpy.utility import check_bounds, scale_samples

[docs]class LHS: @beartype def __init__( self, ndim: int, bounds: Union[np.ndarray, None] = None, seed: Union[int, None] = None, criterion: str = 'random', iteration: Union[int, None] = None) -> None: """Latin hypercube sampling. Parameters ---------- ndim : int Dimension of parameters. bounds : numpy array Upper and lower boundaries of each parameter. Shape :code:`(ndim, 2)`. `bounds[:, 0]` corresponds to lower boundaries of each parameter and `bounds[:, 1]` to upper boundaries of each parameter. seed : int, optional Seed to initialize the pseudo-random number generator. criterion : str, optional Criterion for generating Latin hypercube samples. `'random'` - randomly locate samples in each bin. `'center'` - locate samples in bin centers. `'maximin'` - locate samples in each bin by maximizing their minimum distance. Default is `'random'`. iteration : int, optional Number of iterations if :code:`criterion='maximin'`. """ self.ndim = ndim if bounds is not None: check_bounds(ndim, bounds) self.bounds = bounds self.rng = np.random.default_rng(seed) self.seed = seed if criterion not in ('random', 'center', 'maximin'): raise NotImplementedError( f"criterion {criterion} is not implemented." f" Supported criterion includes: 'random', 'center', 'maximin'") self.criterion = criterion self.iteration = iteration
[docs] def sample(self, nsamples: int) -> np.ndarray: """Draw Latin hypercube samples. Parameters ---------- nsamples : int Number of samples to be drawn. Returns ------- lhs_samples : numpy array Latin hypercube samples. Shape :code:`(nsamples, ndim)`. """ self.nsamples = nsamples # create nsamples bins in the interval [0, 1] bin_edges = np.linspace(0, 1, nsamples+1) self._lower_edges = bin_edges[:nsamples] self._upper_edges = bin_edges[1 : nsamples+1] if self.criterion == 'random': lhs_samples = self._random() elif self.criterion == 'center': lhs_samples = self._center() elif self.criterion == 'maximin': lhs_samples = self._maximin() if self.bounds is not None: lhs_samples = scale_samples(lhs_samples, self.bounds) return lhs_samples
def _random(self) -> np.ndarray: """"Draw samples by randomly picking points in each bin.""" samples = self.rng.uniform(size=(self.nsamples, self.ndim)) lhs_samples = np.zeros_like(samples) for j in range(self.ndim): samples[:,j] = samples[:,j]*(self._upper_edges-self._lower_edges) \ + self._lower_edges lhs_samples[:,j] = self.rng.permutation(samples[:,j]) return lhs_samples def _center(self) -> np.ndarray: """"Draw samples by picking points at the center of each bin.""" self._center_points = (self._lower_edges + self._upper_edges) / 2 lhs_samples = np.zeros((self.nsamples, self.ndim)) for j in range(self.ndim): lhs_samples[:,j] = self.rng.permutation(self._center_points) return lhs_samples def _maximin(self) -> np.ndarray: """" Draw samples by maximizing the minimum distance between samples. """ if self.iteration is None: self.iteration = 100 # max --> min? min_distance = 0 lhs_samples = np.zeros((self.nsamples, self.ndim)) for i in range(self.iteration): temp_samples = self._random() temp_distance = pdist(temp_samples) if np.min(temp_distance) > min_distance: min_distance = np.min(temp_distance) lhs_samples = temp_samples return lhs_samples