diff --git a/doc/Makefile b/doc/Makefile
index 195e17ce0a2958904630a731b51ecfe4204d129c..4611b6aee9cc85babbc81d75af2c5e888534fed7 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -3,7 +3,7 @@
# You can set these variables from the command line.
SPHINXOPTS =
-SPHINXBUILD = sphinx-build
+SPHINXBUILD = python $(shell which sphinx-build)
SPHINXPROJ = pytools
SOURCEDIR = .
BUILDDIR = _build
@@ -17,4 +17,4 @@ help:
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
- @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/pytools/__init__.py b/pytools/__init__.py
index 64a896bd4540cededbd2f4a015356dba8768d7d6..4ee4fc124bee7a347aa9cbc92d6dc85191ba0498 100644
--- a/pytools/__init__.py
+++ b/pytools/__init__.py
@@ -37,7 +37,6 @@ from typing import (
Any, Callable, Dict, Iterable, List, Optional, Set, Tuple, TypeVar)
import builtins
-
from sys import intern
@@ -171,7 +170,8 @@ Hashing
Sampling
--------
-.. autofunction:: sphere_sample
+.. autofunction:: sphere_sample_equidistant
+.. autofunction:: sphere_sample_fibonacci
Type Variables Used
-------------------
@@ -256,7 +256,7 @@ def deprecate_keyword(oldkey: str,
# }}}
-# {{{ math --------------------------------------------------------------------
+# {{{ math
def delta(x, y):
if x == y:
@@ -472,7 +472,7 @@ class FakeList:
return self._Function(index)
-# {{{ dependent dictionary ----------------------------------------------------
+# {{{ dependent dictionary
class DependentDictionary:
def __init__(self, f, start=None):
@@ -1809,7 +1809,7 @@ def word_wrap(text, width, wrap_using="\n"):
# }}}
-# {{{ command line interfaces -------------------------------------------------
+# {{{ command line interfaces
def _exec_arg(arg, execenv):
import os
@@ -2701,11 +2701,12 @@ def unordered_hash(hash_instance, iterable, hash_constructor=None):
# {{{ sphere_sample
-def sphere_sample(npoints_approx, r=1):
- """Generate points regularly distributed on a sphere.
- Based on: https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
- Returns an array of shape (npoints, 3). npoints does not generally equally
- npoints_approx.
+def sphere_sample_equidistant(npoints_approx: int, r: float = 1.0):
+ """Generate points regularly distributed on a sphere
+ based on https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf.
+
+ :returns: an :class:`~numpy.ndarray` of shape ``(3, npoints)``, where
+ ``npoints`` does not generally equal *npoints_approx*.
"""
import numpy as np
@@ -2728,14 +2729,59 @@ def sphere_sample(npoints_approx, r=1):
r * np.cos(theta)])
count += 1
- # Add poles.
+ # add poles
for i in range(3):
for sign in [-1, +1]:
pole = np.zeros(3)
pole[i] = r * sign
points.append(pole)
- return np.array(points)
+ return np.array(points).T.copy()
+
+
+# NOTE: each tuple contains ``(epsilon, max_npoints)``
+_SPHERE_FIBONACCI_OFFSET = (
+ (0.33, 24), (1.33, 177), (3.33, 890),
+ (10, 11000), (27, 39000), (75, 600000), (214, float("inf")),
+ )
+
+
+def sphere_sample_fibonacci(
+ npoints: int, r: float = 1.0, *,
+ optimize: Optional[str] = None):
+ """Generate points on a sphere based on an offset Fibonacci lattice from [2]_.
+
+ .. [2] http://extremelearning.com.au/how-to-evenly-distribute-points-on-a-sphere-more-effectively-than-the-canonical-fibonacci-lattice/
+
+ :param optimize: takes the values: *None* to use the standard Fibonacci
+ lattice, ``"minimum"`` to minimize the nearest neighbor distances in the
+ lattice and ``"average"`` to minimize the average distances in the
+ lattice.
+
+ :returns: an :class:`~numpy.ndarray` of shape ``(3, npoints)``.
+ """ # noqa: E501
+
+ import numpy as np
+ if optimize is None:
+ epsilon = 0.5
+ elif optimize == "minimum":
+ epsilon, _ = next(o for o in _SPHERE_FIBONACCI_OFFSET if npoints < o[1])
+ elif optimize == "average":
+ epsilon = 0.36
+ else:
+ raise ValueError(f"unknown 'optimize' choice: '{optimize}'")
+
+ golden_ratio = (1 + np.sqrt(5)) / 2
+ n = np.arange(npoints)
+
+ phi = 2.0 * np.pi * n / golden_ratio
+ theta = np.arccos(1.0 - 2.0 * (n + epsilon) / (npoints + 2 * epsilon - 1))
+
+ return np.stack([
+ r * np.sin(theta) * np.cos(phi),
+ r * np.sin(theta) * np.sin(phi),
+ r * np.cos(theta)
+ ])
# }}}
diff --git a/test/test_pytools.py b/test/test_pytools.py
index 6c24ecc501c80dde26566a48ad86bb8924a25dc5..17ee1fff8fbb425c028d155accece9348ff44004 100644
--- a/test/test_pytools.py
+++ b/test/test_pytools.py
@@ -555,6 +555,61 @@ def test_unordered_hash():
!= unordered_hash(hashlib.sha256(), lst).digest())
+# {{{ sphere sampling
+
+@pytest.mark.parametrize("sampling", [
+ "equidistant", "fibonacci", "fibonacci_min", "fibonacci_avg",
+ ])
+def test_sphere_sampling(sampling, visualize=False):
+ from functools import partial
+ from pytools import sphere_sample_equidistant, sphere_sample_fibonacci
+
+ npoints = 128
+ radius = 1.5
+
+ if sampling == "equidistant":
+ sampling_func = partial(sphere_sample_equidistant, r=radius)
+ elif sampling == "fibonacci":
+ sampling_func = partial(
+ sphere_sample_fibonacci, r=radius, optimize=None)
+ elif sampling == "fibonacci_min":
+ sampling_func = partial(
+ sphere_sample_fibonacci, r=radius, optimize="minimum")
+ elif sampling == "fibonacci_avg":
+ sampling_func = partial(
+ sphere_sample_fibonacci, r=radius, optimize="average")
+ else:
+ raise ValueError(f"unknown sampling method: '{sampling}'")
+
+ import numpy as np
+ points = sampling_func(npoints)
+ assert np.all(np.linalg.norm(points, axis=0) < radius + 1.0e-15)
+
+ if not visualize:
+ return
+
+ import matplotlib.pyplot as plt
+ fig = plt.figure(figsize=(10, 10), dpi=300)
+ ax = fig.add_subplot(111, projection="3d")
+
+ import matplotlib.tri as mtri
+ theta = np.arctan2(np.sqrt(points[0]**2 + points[1]**2), points[2])
+ phi = np.arctan2(points[1], points[0])
+ triangles = mtri.Triangulation(theta, phi)
+
+ ax.plot_trisurf(points[0], points[1], points[2], triangles=triangles.triangles)
+ ax.set_xlim([-radius, radius])
+ ax.set_ylim([-radius, radius])
+ ax.set_zlim([-radius, radius])
+ ax.margins(0.05, 0.05, 0.05)
+
+ # plt.show()
+ fig.savefig(f"sphere_sampling_{sampling}")
+ plt.close(fig)
+
+# }}}
+
+
if __name__ == "__main__":
if len(sys.argv) > 1:
exec(sys.argv[1])