Repository: scipy/scipy
Pull request: ENH: Geometric quality metrics for random sampling 
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Depending on your changes, you can skip CI operations and save time and energy: 
http://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
Closes #14037

#### What does this implement/fix?
The change adds simple geometric/topological metrics for random sampling.

#### Additional information
The metrics are implemented in a new function rather than added to the existing `scipy.stats.qmc.discrepancy` function as they represent a different type of metrics: distance/geometric-based vs distribution-based. The interpretation of the returned value is also different between the two functions. See the discussion in the comments below (particularly, https://github.com/scipy/scipy/pull/18893#pullrequestreview-1531761774 and https://github.com/scipy/scipy/pull/18893#issuecomment-1650105609) for further details.


Problem statement:
ENH: other quality metrics for random sampling
The new `scipy.stats.qmc` sub-module provide `discrepancy` as a way to measure the quality of a random sample. This measure is a uniformity criterion.

Other metrics can be used to measure the quality such as _mindist_ and _minimum spanning tree_. 

Here is a simple implementation of useful measures:

```python
import numpy as np
from scipy.spatial.distance import cdist


def l1_norm(sample):
    l1 = cdist(sample, sample, 'cityblock')
    return np.min(l1[l1.nonzero()])


# mindist
def l2_norm(sample):
    l2 = cdist(sample, sample)
    return np.min(l2[l2.nonzero()])


def linf_norm(sample):
    linf = cdist(sample, sample, 'chebyshev')
    return np.min(linf[linf.nonzero()])
```

These are more "geometrically/topologically" based whereas the discrepancy is a statistical criterion. We could have a common function with a `method` argument to select one of the metrics.

From a dependency perspective, this should go in `scipy.spatial`. Also, if we add _mst_, `scipy.spatial` already has a dependency to `scipy.sparse`.

Hints from issue discussion:
@tupui Can I take this up too? Are there any discussions needed for deciding which measures are to be included? Related to https://github.com/scipy/scipy/issues/14510
> @tupui Can I take this up too? Are there any discussions needed for deciding which measures are to be included? Related to #14510

Hey @czgdp1807, this one is free for all, but very accessible this is why I flagged it for sprints. So I would prefer if you could find another one to help on (like poisson disk as I wrote in another issue). But OC if you would really like to take something I would not reject it 😄 It should be fun to come here and help out.
Can I work on this as part of the today's sprint?
@kejsitake I am happy to review any work there. Warning though this is quite involved (needs docs, tests, etc.) and would also need to go through the mailing list at some point. So it would probably not make it in the 1.10 release (also because I will only have this week to review and then I am back in Dec. right before the branch for the new release)
@kshitijaggarwal and myself would like to take this up for the SciPy 2023 sprint.
I am happy to review anything here.

@mdhaber what do you think about making a single function for these? Something like `geometric_measures` or `geometric_uniformity`?
I don't object. I can't personally envision usage of the function so I can't be more assertive : )
While working on an implementation, we found that the `minimum_spanning_tree` function does not distinguish between the absence of an edge and an edge with zero weight in its input: https://github.com/scipy/scipy/issues/18892.

A work-in-progress pull request is here for any early feedback: https://github.com/scipy/scipy/pull/18893

We had a few questions:

- Should we add this function in a new module in `scipy.spatial`, keep it in `scipy.stats` or perhaps merge it into `scipy.stats.qmc.discrepancy`?
- We can get two metrics from the minimum spanning tree: the mean distance and the standard deviation, so we could let the user choose which metric they want or return both.

Thank you 😃 I will reply on the PR itself.

Commit message: TST: Tests for geometric_discrepancy

Patch to review:
diff --git a/scipy/stats/_qmc.py b/scipy/stats/_qmc.py
index 647948c05e06..6f680b68b787 100644
--- a/scipy/stats/_qmc.py
+++ b/scipy/stats/_qmc.py
@@ -26,6 +26,7 @@
 
 import scipy.stats as stats
 from scipy._lib._util import rng_integers, _rng_spawn
+from scipy.sparse.csgraph import minimum_spanning_tree
 from scipy.spatial import distance, Voronoi
 from scipy.special import gammainc
 from ._sobol import (
@@ -43,7 +44,7 @@
 )
 
 
-__all__ = ['scale', 'discrepancy', 'update_discrepancy',
+__all__ = ['scale', 'discrepancy', 'geometric_discrepancy', 'update_discrepancy',
            'QMCEngine', 'Sobol', 'Halton', 'LatinHypercube', 'PoissonDisk',
            'MultinomialQMC', 'MultivariateNormalQMC']
 
@@ -168,6 +169,36 @@ def scale(
         return (sample - lower) / (upper - lower)
 
 
+def _ensure_in_unit_hypercube(sample: npt.ArrayLike) -> np.ndarray:
+    """Ensure that sample is a 2D array and is within a unit hypercube
+
+    Parameters
+    ----------
+    sample : array_like (n, d)
+        A 2D array of points.
+
+    Returns
+    -------
+    np.ndarray
+        The array interpretation of the input sample
+
+    Raises
+    ------
+    ValueError
+        If the input is not a 2D array or contains points outside of
+        a unit hypercube.
+    """
+    sample = np.asarray(sample, dtype=np.float64, order="C")
+
+    if not sample.ndim == 2:
+        raise ValueError("Sample is not a 2D array")
+
+    if (sample.max() > 1.) or (sample.min() < 0.):
+        raise ValueError("Sample is not in unit hypercube")
+
+    return sample
+
+
 def discrepancy(
         sample: npt.ArrayLike,
         *,
@@ -195,6 +226,10 @@ def discrepancy(
     discrepancy : float
         Discrepancy.
 
+    See Also
+    --------
+    geometric_discrepancy
+
     Notes
     -----
     The discrepancy is a uniformity criterion used to assess the space filling
@@ -281,14 +316,7 @@ def discrepancy(
     0.008142039609053513
 
     """
-    sample = np.asarray(sample, dtype=np.float64, order="C")
-
-    # Checking that sample is within the hypercube and 2D
-    if not sample.ndim == 2:
-        raise ValueError("Sample is not a 2D array")
-
-    if (sample.max() > 1.) or (sample.min() < 0.):
-        raise ValueError("Sample is not in unit hypercube")
+    sample = _ensure_in_unit_hypercube(sample)
 
     workers = _validate_workers(workers)
 
@@ -306,6 +334,128 @@ def discrepancy(
                          f" {set(methods)!r}")
 
 
+def geometric_discrepancy(
+        sample: npt.ArrayLike,
+        method: Literal["mindist", "mst"] = "mindist",
+        metric: str = "euclidean") -> float:
+    """Discrepancy of a given sample based on its geometric properties.
+
+    Parameters
+    ----------
+    sample : array_like (n, d)
+        The sample to compute the discrepancy from.
+    method : {"mindist", "mst"}, optional
+        The method to use. One of ``mindist`` for minimum distance (default)
+        or ``mst`` for minimum spanning tree.
+    metric : str or callable, optional
+        The distance metric to use. See the documentation
+        for `scipy.spatial.distance.pdist` for the available metrics and
+        the default.
+
+    Returns
+    -------
+    discrepancy : float
+        Discrepancy (higher values correspond to greater sample uniformity).
+
+    See Also
+    --------
+    discrepancy
+
+    Notes
+    -----
+    The discrepancy can serve as a simple measure of quality of a random sample.
+    This measure is based on the geometric properties of the distribution of points
+    in the sample, such as the minimum distance between any pair of points, or
+    the mean edge length in a minimum spanning tree.
+
+    The higher the value is, the better the coverage of the parameter space is.
+    Note that this is different from `scipy.stats.qmc.discrepancy`, where lower
+    values correspond to higher quality of the sample.
+
+    It is possible to calculate two metrics from the minimum spanning tree:
+    the mean edge length and the standard deviation of edges lengths. Using
+    both metrics offers a better picture of uniformity than either metric alone,
+    with higher mean and lower standard deviation being preferable (see [1]_
+    for a brief discussion). This function currently only calculates the mean
+    edge length.
+
+    References
+    ----------
+    .. [1] Franco J. et al. "Minimum Spanning Tree: A new approach to assess the quality
+       of the design of computer experiments." Chemometrics and Intelligent Laboratory
+       Systems, 97 (2), pp. 164-169, 2009.
+
+    Examples
+    --------
+    Calculate the quality of the sample using the minimum euclidean distance
+    (the defaults):
+
+    >>> import numpy as np
+    >>> from scipy.stats import qmc
+    >>> rng = np.random.default_rng(191468432622931918890291693003068437394)
+    >>> sample = qmc.LatinHypercube(d=2, seed=rng).random(50)
+    >>> qmc.geometric_discrepancy(sample)
+    0.03708161435687876
+
+    Calculate the quality using the mean edge length in the minimum
+    spanning tree:
+
+    >>> qmc.geometric_discrepancy(sample, method='mst')
+    0.1105149978798376
+
+    Display the minimum spanning tree and the points with
+    the smallest distance:
+
+    >>> import matplotlib.pyplot as plt
+    >>> from matplotlib.lines import Line2D
+    >>> from scipy.sparse.csgraph import minimum_spanning_tree
+    >>> from scipy.spatial.distance import pdist, squareform
+    >>> dist = pdist(sample)
+    >>> mst = minimum_spanning_tree(squareform(dist))
+    >>> edges = np.where(mst.toarray() > 0)
+    >>> edges = np.asarray(edges).T
+    >>> min_dist = np.min(dist)
+    >>> min_idx = np.argwhere(squareform(dist) == min_dist)[0]
+    >>> fig, ax = plt.subplots(figsize=(10, 5))
+    >>> _ = ax.set(aspect='equal', xlabel=r'$x_1$', ylabel=r'$x_2$',
+    ...            xlim=[0, 1], ylim=[0, 1])
+    >>> for edge in edges:
+    ...     ax.plot(sample[edge, 0], sample[edge, 1], c='k')
+    >>> ax.scatter(sample[:, 0], sample[:, 1])
+    >>> ax.add_patch(plt.Circle(sample[min_idx[0]], min_dist, color='red', fill=False))
+    >>> markers = [
+    ...     Line2D([0], [0], marker='o', lw=0, label='Sample points'),
+    ...     Line2D([0], [0], color='k', label='Minimum spanning tree'),
+    ...     Line2D([0], [0], marker='o', lw=0, markerfacecolor='w', markeredgecolor='r',
+    ...            label='Minimum point-to-point distance'),
+    ... ]
+    >>> ax.legend(handles=markers, loc='center left', bbox_to_anchor=(1, 0.5));
+    >>> plt.show()
+
+    """
+    sample = _ensure_in_unit_hypercube(sample)
+    if sample.shape[0] < 2:
+        raise ValueError("Sample must contain at least two points")
+
+    distances = distance.pdist(sample, metric=metric)  # type: ignore[call-overload]
+
+    if np.any(distances == 0.0):
+        warnings.warn("Sample contains duplicate points.", stacklevel=2)
+
+    if method == "mindist":
+        return np.min(distances[distances.nonzero()])
+    elif method == "mst":
+        fully_connected_graph = distance.squareform(distances)
+        mst = minimum_spanning_tree(fully_connected_graph)
+        distances = mst[mst.nonzero()]
+        # TODO consider returning both the mean and the standard deviation
+        # see [1] for a discussion
+        return np.mean(distances)
+    else:
+        raise ValueError(f"{method!r} is not a valid method. "
+                         f"It must be one of {{'mindist', 'mst'}}")
+
+
 def update_discrepancy(
         x_new: npt.ArrayLike,
         sample: npt.ArrayLike,
diff --git a/scipy/stats/qmc.py b/scipy/stats/qmc.py
index cabb7d5b103d..337f4987ec1d 100644
--- a/scipy/stats/qmc.py
+++ b/scipy/stats/qmc.py
@@ -33,6 +33,7 @@
    :toctree: generated/
 
    discrepancy
+   geometric_discrepancy
    update_discrepancy
    scale
 
diff --git a/scipy/stats/tests/test_qmc.py b/scipy/stats/tests/test_qmc.py
index 2fc1cb0c355a..d351ac45f6f3 100644
--- a/scipy/stats/tests/test_qmc.py
+++ b/scipy/stats/tests/test_qmc.py
@@ -178,6 +178,56 @@ def test_discrepancy_parallel(self, monkeypatch):
         with pytest.raises(ValueError, match="Invalid number of workers..."):
             qmc.discrepancy(sample, workers=-2)
 
+    def test_geometric_discrepancy_errors(self):
+        sample = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+
+        with pytest.raises(ValueError, match=r"Sample is not in unit hypercube"):
+            qmc.geometric_discrepancy(sample)
+
+        with pytest.raises(ValueError, match=r"Sample is not a 2D array"):
+            qmc.geometric_discrepancy([1, 3])
+
+        sample = [[0, 0], [1, 1], [0.5, 0.5]]
+        with pytest.raises(ValueError, match=r"'toto' is not a valid ..."):
+            qmc.geometric_discrepancy(sample, method="toto")
+
+        sample = np.array([[0, 0], [0, 0], [0, 1]])
+        with pytest.raises(UserWarning, match="Sample contains duplicate points."):
+            qmc.geometric_discrepancy(sample)
+
+        sample = np.array([[0.5, 0.5]])
+        with pytest.raises(ValueError, match="Sample must contain at least two points"):
+            qmc.geometric_discrepancy(sample)
+
+    def test_geometric_discrepancy(self):
+        sample = np.array([[0, 0], [1, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), np.sqrt(2))
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), np.sqrt(2))
+
+        sample = np.array([[0, 0], [0, 1], [0.5, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), 0.5)
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), 0.75)
+
+        sample = np.array([[0, 0], [0.25, 0.25], [1, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), np.sqrt(2) / 4)
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), np.sqrt(2) / 2)
+        assert_allclose(qmc.geometric_discrepancy(sample, metric="chebyshev"), 0.25)
+        assert_allclose(
+            qmc.geometric_discrepancy(sample, method="mst", metric="chebyshev"), 0.5
+        )
+
+        rng = np.random.default_rng(191468432622931918890291693003068437394)
+        sample = qmc.LatinHypercube(d=3, seed=rng).random(50)
+        assert_allclose(qmc.geometric_discrepancy(sample), 0.05106012076093356)
+        assert_allclose(
+            qmc.geometric_discrepancy(sample, method='mst'), 0.19704396643366182
+        )
+
+    @pytest.mark.xfail(reason="minimum_spanning_tree ignores zero distances (#18892)")
+    def test_geometric_discrepancy_mst_with_zero_distances(self):
+        sample = np.array([[0, 0], [0, 0], [0, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample, method='mst'), 0.5)
+
     def test_update_discrepancy(self):
         # From Fang et al. Design and modeling for computer experiments, 2006
         space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])


Code review comments:
diff_hunk: @@ -178,6 +178,56 @@ def test_discrepancy_parallel(self, monkeypatch):
         with pytest.raises(ValueError, match="Invalid number of workers..."):
             qmc.discrepancy(sample, workers=-2)
 
+    def test_geometric_discrepancy_errors(self):
+        sample = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
+
+        with pytest.raises(ValueError, match=r"Sample is not in unit hypercube"):
+            qmc.geometric_discrepancy(sample)
+
+        with pytest.raises(ValueError, match=r"Sample is not a 2D array"):
+            qmc.geometric_discrepancy([1, 3])
+
+        sample = [[0, 0], [1, 1], [0.5, 0.5]]
+        with pytest.raises(ValueError, match=r"'toto' is not a valid ..."):
+            qmc.geometric_discrepancy(sample, method="toto")
+
+        sample = np.array([[0, 0], [0, 0], [0, 1]])
+        with pytest.raises(UserWarning, match="Sample contains duplicate points."):
+            qmc.geometric_discrepancy(sample)
+
+        sample = np.array([[0.5, 0.5]])
+        with pytest.raises(ValueError, match="Sample must contain at least two points"):
+            qmc.geometric_discrepancy(sample)
+
+    def test_geometric_discrepancy(self):
+        sample = np.array([[0, 0], [1, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), np.sqrt(2))
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), np.sqrt(2))
+
+        sample = np.array([[0, 0], [0, 1], [0.5, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), 0.5)
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), 0.75)
+
+        sample = np.array([[0, 0], [0.25, 0.25], [1, 1]])
+        assert_allclose(qmc.geometric_discrepancy(sample), np.sqrt(2) / 4)
+        assert_allclose(qmc.geometric_discrepancy(sample, method="mst"), np.sqrt(2) / 2)
+        assert_allclose(qmc.geometric_discrepancy(sample, metric="chebyshev"), 0.25)
+        assert_allclose(
+            qmc.geometric_discrepancy(sample, method="mst", metric="chebyshev"), 0.5
+        )
+
+        rng = np.random.default_rng(191468432622931918890291693003068437394)
+        sample = qmc.LatinHypercube(d=3, seed=rng).random(50)
+        assert_allclose(qmc.geometric_discrepancy(sample), 0.05106012076093356)
+        assert_allclose(
+            qmc.geometric_discrepancy(sample, method='mst'), 0.19704396643366182
+        )
+
+    @pytest.mark.xfail(reason="minimum_spanning_tree ignores zero distances (#18892)")
+    def test_geometric_discrepancy_mst_with_zero_distances(self):
path: scipy/stats/tests/test_qmc.py
text: @user2:
I like that, this is a test which should pass once the issue is resolve correct? If that's the case, then we can add the strict parameter which would make the test fail once this pass.

```suggestion
    @user1.mark.xfail(
        reason="minimum_spanning_tree ignores zero distances (#18892)",
        strict=True
    )
    def test_geometric_discrepancy_mst_with_zero_distances(self):
```

@author:
Done