# Source code for psdr.geometry.fill

from __future__ import print_function, division

import numpy as np
from scipy.spatial.distance import cdist

from .vertex import voronoi_vertex_sample, voronoi_vertex

def fill_distance(domain, Xhat, L = None):
r""" Compute the exact fill distance by constructing the bounded Voronoi vertices

Parameters
----------
domain: Domain
Domain on which to compute the dispersion
Xhat: array-like (?, m)
Existing samples on the domain
L: array-like (?, m) optional
Matrix defining the distance metric on the domain

Returns
-------
d: float
Fill distance
"""

V = voronoi_vertex(domain, Xhat, L = L)
D = psdr.cdist(X, V, L = L)
return np.max(np.min(D, axis= 0))

[docs]def fill_distance_estimate(domain, Xhat, L = None, Nsamp = int(1e3), X0 = None ):
r""" Estimate the fill distance of the points Xhat in the domain

The *fill distance* (Def. 1.4 of [Wen04]_) or *dispersion* [LC05]_
is the furthest distance between any point :math:\mathbf{x} \in \mathcal{D}
and a set of points
:math:\lbrace \widehat{\mathbf{x}}_j \rbrace_{j=1}^m \subset \mathcal{D}:

.. math::

\sup_{\mathbf{x} \in \mathcal{D}} \min_{j=1,\ldots,M} \|\mathbf{L}(\mathbf{x} - \widehat{\mathbf{x}}_j)\|_2.

Similar to :meth:psdr.seq_maximin_sample, this uses :meth:psdr.voronoi_vertex_sample to find
a subset of local maximizers and returns the best of these.

Parameters
----------
domain: Domain
Domain on which to compute the dispersion
Xhat: array-like (?, m)
Existing samples on the domain
L: array-like (?, m) optional
Matrix defining the distance metric on the domain
Nsamp: int, default 1e4
Number of samples to use for vertex sampling
X0: array-like (?, m)
Samples from the domain to use in :meth:psdr.voronoi_vertex_sample

Returns
-------
d: float
Fill distance lower bound

References
----------
..  [Wen04] Scattered Data Approximation. Holger Wendland.
Cambridge University Press, 2004.
https://doi.org/10.1017/CBO9780511617539
.. [LC05] Iteratively Locating Voronoi Vertices for Dispersion Estimation
Stephen R. Lindemann and Peng Cheng
Proceedings of the 2005 Interational Conference on Robotics and Automation
"""
if L is None:
L = np.eye(len(domain))

if X0 is None:
from ..sample import initial_sample
X0 = initial_sample(domain, L, Nsamp = Nsamp)

# Since we only care about distance in L, we can terminate early if L is rank-deficient
# and hence we turn off the randomization
Xcan = voronoi_vertex_sample(domain, Xhat, X0, L = L, randomize = False)

# Euclidean distance
D = cdist(L.dot(Xcan.T).T, L.dot(Xhat.T).T)

d = np.min(D, axis = 1)
return float(np.max(d))