Repository: scipy/scipy
Pull request: ENH: Poisson Disk sampling for QMC
#### Reference issue
Adding Poisson Disk sampling to QMC #13912

#### What does this implement/fix?

Poisson disk sampling in arbitrary dimensions using Bridson's algorithm, implemented using ``numpy`` and ``scipy``.
Generates so-called "blue noise" that prevents clustering by ensuring each two points are at least ``radius`` apart.
https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

<details>

```python
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from scipy.stats import qmc


rng = np.random.default_rng()
radius = 0.2
engine = qmc.PoissonDisk(d=2, radius=0.2, seed=rng)
sample = engine.random(20)

fig, ax = plt.subplots()
_ = ax.scatter(sample[:, 0], sample[:, 1])
circles = [plt.Circle((xi, yi), radius=radius/2, fill=False)
           for xi, yi in sample]
collection = PatchCollection(circles, match_original=True)
ax.add_collection(collection)
_ = ax.set(aspect='equal', xlabel=r'$x_1$', ylabel=r'$x_2$',
           xlim=[0, 1], ylim=[0, 1])
plt.show()
```
</details>

![poissonDisk](https://user-images.githubusercontent.com/23188539/169402582-5331019f-9093-4456-8a40-94c900de6c9c.png)



Problem statement:
Adding Poisson Disc sampling to QMC
Hi @tupui,

I have looked at the class definition for QMC engines, and I'm having some trouble adapting Bridson's algorithm for PD to these definitions.

As far as I see, I should be able to produce an arbitrary number of samples when `random (n)` is called. There are two issues there due to the nature of the distribution:

1. Since every new sample excludes a sphere around itself from any further sampling, there is a finite number of samples that can be drawn from a region of space for a given minimum distance between points _r_. What should the function do if it can produce no more samples?
2. Bridson's algorithm exploits the properties of a maximum PD distribution to produce new points efficiently. Specifically, it draws new points only in an annulus _[r, 2r]_ away from an existing point. This means, that if we draw PD points sequentially, we will densely fill one region of space first, instead of drawing uniformly in the [0; 1) interval. One way that I see to circumvent this is to generate the entire distribution on `__init__`, and then to return random samples from this distribution each time `random(n)` is called. Is this approach acceptable for the library?

Hints from issue discussion:
Thanks for tackling this! Poisson disk sampling is quite useful and would definitely be valuable.

I am not certain yet about the final interface. For 1. we could raise a warning. And for 2. We could also do like with `LatinHypercube` where calling again `random` would produce a new sampling (uncorrelated from the first one). We can think later about fancy strategies to add points (I've seen things like varying radius, etc. Same for LHS). But for now I would just propose a clean implementation of the base sampler. This is already valuable, especially if you have something in n-dimensions.
For reference, xref https://github.com/scikit-image/scikit-image/issues/2380
FYI, feel free to create a PR in draft mode if you need help/advice. You don't need to have something totally finished to propose it. I/we can help you to figure out things 😃
Thanks! I've gone ahead and created a draft PR. Right now it's very basic, e.g. I don't use `centered` or `seed` in any way, and each call of `random` produces a new uncorrelated sample, with the issues I've outlined above

Commit message: MAINT: refactor code.

[skip ci]

Patch to review:
diff --git a/scipy/stats/_qmc.py b/scipy/stats/_qmc.py
index fcc11dd92347..41d1d340711a 100644
--- a/scipy/stats/_qmc.py
+++ b/scipy/stats/_qmc.py
@@ -30,6 +30,7 @@
 import scipy.stats as stats
 from scipy._lib._util import rng_integers
 from scipy.spatial import distance, Voronoi
+from scipy.special.cython_special import gammainc
 from ._sobol import (
     _initialize_v, _cscramble, _fill_p_cumulative, _draw, _fast_forward,
     _categorize, _MAXDIM
@@ -1595,6 +1596,183 @@ def fast_forward(self, n: IntNumber) -> Sobol:
         return self
 
 
+class PoissonDisc(QMCEngine):
+    """Poisson disk sampling.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    hypersphere : {"volume", "surface"}, optional
+        Sampling strategy to generate potential candidates to be added in the
+        final sample. Default is "volume".
+
+        * ``volume``: original Bridson algorithm as described in [1]_.
+          new candidates are sampled within the hypersphere.
+        * ``surface``: only sample the surface of the hypersphere.
+
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``Generator`` instance then that instance is
+        used.
+
+    Notes
+    -----
+
+    Taken from [2]_ by P. Zun, written consent given on 31.03.2021
+    by the original author, Shamis, for free use in SciPy under
+    the 3-clause BSD.
+
+    References
+    ----------
+    .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
+       Dimensions." SIGGRAPH, 2007.
+    .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`_.
+
+    """
+
+    def __init__(
+        self,
+        d: IntNumber,
+        hypersphere: Literal["volume", "surface"] = "volume",
+        seed: SeedType = None
+    ) -> None:
+        super().__init__(d=d, seed=seed)
+
+        hypersphere_sample = {
+            "volume": self._hypersphere_volume_sample,
+            "surface": self._hypersphere_surface_sample
+        }
+
+        try:
+            self.hypersphere_method = hypersphere_sample[hypersphere]
+        except KeyError as exc:
+            message = (
+                f"{hypersphere!r} is not a valid hypersphere sampling"
+                f" method. It must be one of {set(hypersphere_sample)!r}")
+            raise ValueError(message) from exc
+
+        # size of the sphere from which the samples are drawn relative to the
+        # size of a disc (radius)
+        # for the surface sampler, all new points are almost exactly 1 radius
+        # away from at least one existing sample
+        # eps to avoid rejection
+        self.radius_factor = 2 if hypersphere == "volume" else 1.001
+        self.radius = 0.05
+        self.squared_radius = self.radius**2
+
+        # sample to generate per iteration in the hypersphere around center
+        self.k = 30
+
+        self.cell_size = self.radius / np.sqrt(self.d)
+        self.grid_size = (np.ceil(np.ones(self.d)/self.cell_size)).astype(int)
+
+    def random(self, n: IntNumber = 1) -> np.ndarray:
+        """Draw `n` in the interval ``[0, 1]``.
+
+        Parameters
+        ----------
+        n : int, optional
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            QMC sample.
+
+        """
+
+        def in_limits(sample: np.ndarray) -> bool:
+            return (sample.max() <= 1.) and (sample.min() >= 0.)
+
+        def in_neighborhood(sample: np.ndarray, n: int = 2) -> bool:
+            """
+            Check if there are samples closer than ``squared_radius`` to the
+            candidate `sample`.
+            """
+            indices = (sample / self.cell_size).astype(int)
+            ind_min = np.maximum(indices - n, np.zeros(self.d, dtype=int))
+            ind_max = np.minimum(indices + n + 1, self.grid_size)
+
+            # Check if the center cell is empty
+            if not np.isnan(sample[tuple(indices)][0]):
+                return True
+
+            a = [[slice(ind_min[i], ind_max[i])] for i in range(self.d)]
+
+            with np.errstate(invalid='ignore'):
+                if np.any(
+                    np.sum(
+                        np.square(sample - sample[tuple(a)]), axis=self.d
+                    ) < self.squared_radius
+                ):
+                    return True
+
+            return False
+
+        def add_sample(sample: np.ndarray) -> None:
+            sample_pool.append(sample)
+            indices = (sample / self.cell_size).astype(int)
+            sample[tuple(indices)] = sample
+
+        # Positions of cells
+        # n-dim value for each grid cell
+        sample = np.empty(np.append(self.grid_size, self.d), dtype=np.float32)
+        # Initialise empty cells with NaNs
+        sample.fill(np.nan)
+
+        # the pool is being initialized with a single random sample
+        sample_pool = []
+        add_sample(self.rng.random((1, self.d)))
+        num_drawn = 1
+
+        # exhaust sample pool to have up to n sample
+        while len(sample_pool) and num_drawn < n:
+            # select a sample from the available pool
+            idx_center = self.rng.integers(len(sample_pool))
+            center = sample_pool[idx_center]
+            del sample_pool[idx_center]
+
+            # generate candidates around the center sample
+            candidates = self.hypersphere_method(
+                center, self.radius * self.radius_factor, self.k
+            )
+
+            # keep candidates that satisfy some conditions
+            for candidate in candidates:
+                if in_limits(candidate) and not in_neighborhood(candidate):
+                    add_sample(candidate)
+                    num_drawn = num_drawn + 1
+
+        # filter out NaN from grid
+        sample = sample[~np.isnan(sample).any(axis=self.d)]
+        self.num_generated = num_drawn
+        return sample
+
+    def _hypersphere_volume_sample(
+        self, center: np.ndarray, radius: float, k: int = 1
+    ) -> np.ndarray:
+        """Uniform sampling within hypersphere."""
+        # should remove samples within r/2
+        x = self.rng.standard_normal(size=(k, self.d))
+        ssq = np.sum(x**2, axis=1)
+        fr = radius * gammainc(self.d/2, ssq/2)**(1/self.d) / np.sqrt(ssq)
+        fr_tiled = np.tile(fr.reshape(k, 1), (1, self.d))
+        p = center + np.multiply(x, fr_tiled)
+        return p
+
+    def _hypersphere_surface_sample(
+        self, center: np.ndarray, radius: float, k: int = 1
+    ) -> np.ndarray:
+        """Uniform sampling on the hypersphere's surface."""
+        vec = self.rng.standard_normal(size=(k, self.d))
+        vec /= np.linalg.norm(vec, axis=1)[:, None]
+        p = center + np.multiply(vec, radius)
+        return p
+
+
 class MultivariateNormalQMC:
     r"""QMC sampling from a multivariate Normal :math:`N(\mu, \Sigma)`.
 


Code review comments:
diff_hunk: @@ -1597,117 +1597,180 @@ def fast_forward(self, n: IntNumber) -> Sobol:
 
 
 class PoissonDisc(QMCEngine):
+    """Poisson disk sampling.
 
-    # Uniform sampling in a hypersphere
-    # Based on Matlab implementation by Roger Stafford
-    # Can be optimized for Bridson algorithm by excluding all points within the r/2 sphere
-    def __hypersphere_volume_sample(center,radius,k=1):
-        ndim = center.size
-        x = np.random.normal(size=(k, ndim))
-        ssq = np.sum(x**2,axis=1)
-        fr = radius*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)
-        frtiled = np.tile(fr.reshape(k,1),(1,ndim))
-        p = center + np.multiply(x,frtiled)
-        return p
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    hypersphere : {"volume", "surface"}, optional
+        Sampling strategy to generate potential candidates to be added in the
+        final sample. Default is "volume".
 
-    # Uniform sampling on the sphere's surface
-    def __hypersphere_surface_sample(center,radius,k=1):
-        ndim = center.size
-        vec = np.random.standard_normal(size=(k, ndim))
-        vec /= np.linalg.norm(vec, axis=1)[:,None]
-        p = center + np.multiply(vec, radius)
-        return p
+        * ``volume``: original Bridson algorithm as described in [1]_.
+          new candidates are sampled within the hypersphere.
+        * ``surface``: only sample the surface of the hypersphere.
 
-    def __init__(self, d, centered=False, seed=None):
-        super().__init__(d=d, seed=seed)
-        self.centered = centered
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``Generator`` instance then that instance is
+        used.
 
-    # def Bridson_sampling(dims=np.array([1.0,1.0]), radius=0.05, k=30, hypersphere_sample=hypersphere_volume_sample):
-    def random(self, n=1):
-        # References: Fast Poisson Disk Sampling in Arbitrary Dimensions
-        #             Robert Bridson, SIGGRAPH, 2007
+    Notes
+    -----
 
-        dims=np.ones(self.d)
-        ndim=dims.size
-        numdrawn=0
+    Taken from [2]_ by P. Zun, written consent given on 31.03.2021
+    by the original author, Shamis, for free use in SciPy under
+    the 3-clause BSD.
+
+    References
+    ----------
+    .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
+       Dimensions." SIGGRAPH, 2007.
+    .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`_.
 
-        # these should be parameters set e.g. when the class is declared
-        radius=0.05
-        k=30
+    """
+
+    def __init__(
+        self,
+        d: IntNumber,
+        hypersphere: Literal["volume", "surface"] = "volume",
+        seed: SeedType = None
+    ) -> None:
+        super().__init__(d=d, seed=seed)
 
-        # __hypersphere_volume_sample is the default: classic Bridson algorithm as described here:
-        # https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
-        # __hypersphere_surface_sample is an alternative sampler, results in denser points:
-        # based on the method proposed here by Martin Roberts: http://extremelearning.com.au/an-improved-version-of-bridsons-algorithm-n-for-poisson-disc-sampling/
-        hypersphere_sample=__hypersphere_volume_sample
+        hypersphere_sample = {
+            "volume": self._hypersphere_volume_sample,
+            "surface": self._hypersphere_surface_sample
+        }
 
-        # size of the sphere from which the samples are drawn relative to the size of a disc (radius)
-        sample_factor = 2
-        if hypersphere_sample == __hypersphere_volume_sample:
-            sample_factor = 2
+        try:
+            self.hypersphere_method = hypersphere_sample[hypersphere]
+        except KeyError as exc:
+            message = (
+                f"{hypersphere!r} is not a valid hypersphere sampling"
+                f" method. It must be one of {set(hypersphere_sample)!r}")
+            raise ValueError(message) from exc
 
-        # for the surface sampler, all new points are almost exactly 1 radius away from at least one existing sample
+        # size of the sphere from which the samples are drawn relative to the
+        # size of a disc (radius)
+        # for the surface sampler, all new points are almost exactly 1 radius
+        # away from at least one existing sample
         # eps to avoid rejection
-        if hypersphere_sample == __hypersphere_surface_sample:
-            eps = 0.001
-            sample_factor = 1 + eps
+        self.radius_factor = 2 if hypersphere == "volume" else 1.001
+        self.radius = 0.05
+        self.squared_radius = self.radius**2
 
-        def in_limits(p):
-            return np.all(np.zeros(ndim) <= p) and np.all(p < dims)
+        # sample to generate per iteration in the hypersphere around center
+        self.k = 30
 
-        # Check if there are samples closer than "squared_radius" to the candidate "p"
-        def in_neighborhood(p, n=2):
-            indices = (p / cellsize).astype(int)
-            indmin = np.maximum(indices - n, np.zeros(ndim, dtype=int))
-            indmax = np.minimum(indices + n + 1, gridsize)
+        self.cell_size = self.radius / np.sqrt(self.d)
+        self.grid_size = (np.ceil(np.ones(self.d)/self.cell_size)).astype(int)
+
+    def random(self, n: IntNumber = 1) -> np.ndarray:
+        """Draw `n` in the interval ``[0, 1]``.
+
+        Parameters
+        ----------
+        n : int, optional
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            QMC sample.
+
+        """
+
+        def in_limits(sample: np.ndarray) -> bool:
+            return (sample.max() <= 1.) and (sample.min() >= 0.)
+
+        def in_neighborhood(sample: np.ndarray, n: int = 2) -> bool:
+            """
+            Check if there are samples closer than ``squared_radius`` to the
+            candidate `sample`.
+            """
+            indices = (sample / self.cell_size).astype(int)
+            ind_min = np.maximum(indices - n, np.zeros(self.d, dtype=int))
+            ind_max = np.minimum(indices + n + 1, self.grid_size)
 
             # Check if the center cell is empty
-            if not np.isnan(P[tuple(indices)][0]):
+            if not np.isnan(sample[tuple(indices)][0]):
path: scipy/stats/_qmc.py
text: @author:
This shouldn't be renamed, I think? P large is the array with positions of cells, and p small is the new sample (which you renamed). The same happens in a few other places.

@author:
Ah, sorry. I see this is already renamed in the newest commit.

diff_hunk: @@ -1595,6 +1596,183 @@ def fast_forward(self, n: IntNumber) -> Sobol:
         return self
 
 
+class PoissonDisc(QMCEngine):
+    """Poisson disk sampling.
+
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    hypersphere : {"volume", "surface"}, optional
+        Sampling strategy to generate potential candidates to be added in the
+        final sample. Default is "volume".
+
+        * ``volume``: original Bridson algorithm as described in [1]_.
+          new candidates are sampled within the hypersphere.
+        * ``surface``: only sample the surface of the hypersphere.
+
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``Generator`` instance then that instance is
+        used.
+
+    Notes
+    -----
+
+    Taken from [2]_ by P. Zun, written consent given on 31.03.2021
+    by the original author, Shamis, for free use in SciPy under
+    the 3-clause BSD.
+
+    References
+    ----------
+    .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
+       Dimensions." SIGGRAPH, 2007.
+    .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`_.
+
+    """
+
+    def __init__(
+        self,
+        d: IntNumber,
+        hypersphere: Literal["volume", "surface"] = "volume",
+        seed: SeedType = None
+    ) -> None:
+        super().__init__(d=d, seed=seed)
+
+        hypersphere_sample = {
+            "volume": self._hypersphere_volume_sample,
+            "surface": self._hypersphere_surface_sample
+        }
+
+        try:
+            self.hypersphere_method = hypersphere_sample[hypersphere]
+        except KeyError as exc:
+            message = (
+                f"{hypersphere!r} is not a valid hypersphere sampling"
+                f" method. It must be one of {set(hypersphere_sample)!r}")
+            raise ValueError(message) from exc
+
+        # size of the sphere from which the samples are drawn relative to the
+        # size of a disc (radius)
+        # for the surface sampler, all new points are almost exactly 1 radius
+        # away from at least one existing sample
+        # eps to avoid rejection
+        self.radius_factor = 2 if hypersphere == "volume" else 1.001
+        self.radius = 0.05
+        self.squared_radius = self.radius**2
+
+        # sample to generate per iteration in the hypersphere around center
+        self.k = 30
+
+        self.cell_size = self.radius / np.sqrt(self.d)
+        self.grid_size = (np.ceil(np.ones(self.d)/self.cell_size)).astype(int)
+
+    def random(self, n: IntNumber = 1) -> np.ndarray:
+        """Draw `n` in the interval ``[0, 1]``.
+
+        Parameters
+        ----------
+        n : int, optional
+            Number of samples to generate in the parameter space. Default is 1.
+
+        Returns
+        -------
+        sample : array_like (n, d)
+            QMC sample.
+
+        """
+
+        def in_limits(sample: np.ndarray) -> bool:
+            return (sample.max() <= 1.) and (sample.min() >= 0.)
path: scipy/stats/_qmc.py
text: @user1:
We need this in a few other places in this file already. I would like to create an external function for that.

I will do another PR.

diff_hunk: @@ -1597,117 +1597,180 @@ def fast_forward(self, n: IntNumber) -> Sobol:
 
 
 class PoissonDisc(QMCEngine):
+    """Poisson disk sampling.
 
-    # Uniform sampling in a hypersphere
-    # Based on Matlab implementation by Roger Stafford
-    # Can be optimized for Bridson algorithm by excluding all points within the r/2 sphere
-    def __hypersphere_volume_sample(center,radius,k=1):
-        ndim = center.size
-        x = np.random.normal(size=(k, ndim))
-        ssq = np.sum(x**2,axis=1)
-        fr = radius*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)
-        frtiled = np.tile(fr.reshape(k,1),(1,ndim))
-        p = center + np.multiply(x,frtiled)
-        return p
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    hypersphere : {"volume", "surface"}, optional
+        Sampling strategy to generate potential candidates to be added in the
+        final sample. Default is "volume".
 
-    # Uniform sampling on the sphere's surface
-    def __hypersphere_surface_sample(center,radius,k=1):
-        ndim = center.size
-        vec = np.random.standard_normal(size=(k, ndim))
-        vec /= np.linalg.norm(vec, axis=1)[:,None]
-        p = center + np.multiply(vec, radius)
-        return p
+        * ``volume``: original Bridson algorithm as described in [1]_.
+          new candidates are sampled within the hypersphere.
+        * ``surface``: only sample the surface of the hypersphere.
 
-    def __init__(self, d, centered=False, seed=None):
-        super().__init__(d=d, seed=seed)
-        self.centered = centered
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``Generator`` instance then that instance is
+        used.
 
-    # def Bridson_sampling(dims=np.array([1.0,1.0]), radius=0.05, k=30, hypersphere_sample=hypersphere_volume_sample):
-    def random(self, n=1):
-        # References: Fast Poisson Disk Sampling in Arbitrary Dimensions
-        #             Robert Bridson, SIGGRAPH, 2007
+    Notes
+    -----
 
-        dims=np.ones(self.d)
-        ndim=dims.size
-        numdrawn=0
+    Taken from [2]_ by P. Zun, written consent given on 31.03.2021
+    by the original author, Shamis, for free use in SciPy under
+    the 3-clause BSD.
+
+    References
+    ----------
+    .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
+       Dimensions." SIGGRAPH, 2007.
+    .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`_.
 
-        # these should be parameters set e.g. when the class is declared
-        radius=0.05
-        k=30
+    """
+
+    def __init__(
+        self,
+        d: IntNumber,
+        hypersphere: Literal["volume", "surface"] = "volume",
+        seed: SeedType = None
+    ) -> None:
+        super().__init__(d=d, seed=seed)
 
-        # __hypersphere_volume_sample is the default: classic Bridson algorithm as described here:
-        # https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
-        # __hypersphere_surface_sample is an alternative sampler, results in denser points:
-        # based on the method proposed here by Martin Roberts: http://extremelearning.com.au/an-improved-version-of-bridsons-algorithm-n-for-poisson-disc-sampling/
-        hypersphere_sample=__hypersphere_volume_sample
+        hypersphere_sample = {
+            "volume": self._hypersphere_volume_sample,
+            "surface": self._hypersphere_surface_sample
+        }
 
-        # size of the sphere from which the samples are drawn relative to the size of a disc (radius)
-        sample_factor = 2
-        if hypersphere_sample == __hypersphere_volume_sample:
-            sample_factor = 2
+        try:
+            self.hypersphere_method = hypersphere_sample[hypersphere]
+        except KeyError as exc:
+            message = (
+                f"{hypersphere!r} is not a valid hypersphere sampling"
+                f" method. It must be one of {set(hypersphere_sample)!r}")
+            raise ValueError(message) from exc
 
-        # for the surface sampler, all new points are almost exactly 1 radius away from at least one existing sample
+        # size of the sphere from which the samples are drawn relative to the
+        # size of a disc (radius)
+        # for the surface sampler, all new points are almost exactly 1 radius
+        # away from at least one existing sample
         # eps to avoid rejection
-        if hypersphere_sample == __hypersphere_surface_sample:
-            eps = 0.001
-            sample_factor = 1 + eps
+        self.radius_factor = 2 if hypersphere == "volume" else 1.001
+        self.radius = 0.05
path: scipy/stats/_qmc.py
text: @author:
This should be a parameter for the sampling, it's the minimum distance between the samples

@user1:
Fixed in a later commit

diff_hunk: @@ -1597,117 +1597,180 @@ def fast_forward(self, n: IntNumber) -> Sobol:
 
 
 class PoissonDisc(QMCEngine):
+    """Poisson disk sampling.
 
-    # Uniform sampling in a hypersphere
-    # Based on Matlab implementation by Roger Stafford
-    # Can be optimized for Bridson algorithm by excluding all points within the r/2 sphere
-    def __hypersphere_volume_sample(center,radius,k=1):
-        ndim = center.size
-        x = np.random.normal(size=(k, ndim))
-        ssq = np.sum(x**2,axis=1)
-        fr = radius*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)
-        frtiled = np.tile(fr.reshape(k,1),(1,ndim))
-        p = center + np.multiply(x,frtiled)
-        return p
+    Parameters
+    ----------
+    d : int
+        Dimension of the parameter space.
+    hypersphere : {"volume", "surface"}, optional
+        Sampling strategy to generate potential candidates to be added in the
+        final sample. Default is "volume".
 
-    # Uniform sampling on the sphere's surface
-    def __hypersphere_surface_sample(center,radius,k=1):
-        ndim = center.size
-        vec = np.random.standard_normal(size=(k, ndim))
-        vec /= np.linalg.norm(vec, axis=1)[:,None]
-        p = center + np.multiply(vec, radius)
-        return p
+        * ``volume``: original Bridson algorithm as described in [1]_.
+          new candidates are sampled within the hypersphere.
+        * ``surface``: only sample the surface of the hypersphere.
 
-    def __init__(self, d, centered=False, seed=None):
-        super().__init__(d=d, seed=seed)
-        self.centered = centered
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used,
+        seeded with `seed`.
+        If `seed` is already a ``Generator`` instance then that instance is
+        used.
 
-    # def Bridson_sampling(dims=np.array([1.0,1.0]), radius=0.05, k=30, hypersphere_sample=hypersphere_volume_sample):
-    def random(self, n=1):
-        # References: Fast Poisson Disk Sampling in Arbitrary Dimensions
-        #             Robert Bridson, SIGGRAPH, 2007
+    Notes
+    -----
 
-        dims=np.ones(self.d)
-        ndim=dims.size
-        numdrawn=0
+    Taken from [2]_ by P. Zun, written consent given on 31.03.2021
+    by the original author, Shamis, for free use in SciPy under
+    the 3-clause BSD.
+
+    References
+    ----------
+    .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
+       Dimensions." SIGGRAPH, 2007.
+    .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`_.
 
-        # these should be parameters set e.g. when the class is declared
-        radius=0.05
-        k=30
+    """
+
+    def __init__(
+        self,
+        d: IntNumber,
+        hypersphere: Literal["volume", "surface"] = "volume",
+        seed: SeedType = None
+    ) -> None:
+        super().__init__(d=d, seed=seed)
 
-        # __hypersphere_volume_sample is the default: classic Bridson algorithm as described here:
-        # https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
-        # __hypersphere_surface_sample is an alternative sampler, results in denser points:
-        # based on the method proposed here by Martin Roberts: http://extremelearning.com.au/an-improved-version-of-bridsons-algorithm-n-for-poisson-disc-sampling/
-        hypersphere_sample=__hypersphere_volume_sample
+        hypersphere_sample = {
+            "volume": self._hypersphere_volume_sample,
+            "surface": self._hypersphere_surface_sample
+        }
 
-        # size of the sphere from which the samples are drawn relative to the size of a disc (radius)
-        sample_factor = 2
-        if hypersphere_sample == __hypersphere_volume_sample:
-            sample_factor = 2
+        try:
+            self.hypersphere_method = hypersphere_sample[hypersphere]
+        except KeyError as exc:
+            message = (
+                f"{hypersphere!r} is not a valid hypersphere sampling"
+                f" method. It must be one of {set(hypersphere_sample)!r}")
+            raise ValueError(message) from exc
 
-        # for the surface sampler, all new points are almost exactly 1 radius away from at least one existing sample
+        # size of the sphere from which the samples are drawn relative to the
+        # size of a disc (radius)
+        # for the surface sampler, all new points are almost exactly 1 radius
+        # away from at least one existing sample
         # eps to avoid rejection
-        if hypersphere_sample == __hypersphere_surface_sample:
-            eps = 0.001
-            sample_factor = 1 + eps
+        self.radius_factor = 2 if hypersphere == "volume" else 1.001
+        self.radius = 0.05
+        self.squared_radius = self.radius**2
 
-        def in_limits(p):
-            return np.all(np.zeros(ndim) <= p) and np.all(p < dims)
+        # sample to generate per iteration in the hypersphere around center
+        self.k = 30
path: scipy/stats/_qmc.py
text: @author:
```k``` also should be tunable, I think. 30 is a good guess for most applications (at least in lower dimensions), so it can be a default value. However, in some cases in high dimensions, drawing more candidate samples may be required to have a good distribution. 

@user1:
Ok I can add a parameter. Do you have a better name than `k`? What about `candidates`?

@author:
Yes, ``candidates`` is a good name I think.