Source code for bask.init
import numpy as np
from scipy.optimize import minimize
from sklearn.utils import check_random_state
__all__ = ["sb_sequence", "r2_sequence"]
def _sb_functional(x, X):
if len(X.shape) != 3:
X = X[None, :]
if len(x.shape) == 2:
x = x[:, None, :]
elif len(x.shape) == 1:
x = x[None, None, :]
diff = np.abs(x[..., :] - X)
try:
with np.errstate(divide="raise", invalid="ignore"):
result = np.sum(
np.prod(1 - np.log(2 * np.sin(np.pi * diff)), axis=-1), axis=-1
)
except FloatingPointError:
result = np.inf
return result
def sb_sequence(n, d, existing_points=None, random_state=None, restarts=20):
"""Compute a d-dimensional Steinerberger sequence with up to n points in [0, 1]^d.
If `existing_points` are provided, the remaining points will be generated by the
Steinerberger sequence to fill up the rest of the space.
This function implements the greedy algorithm proposed by Steinerberger 2019 [1]_.
Parameters
----------
n : int
Number of total points to return. If `existing_points` is not None, the sequence
will add points until n is reached.
d : int
Number of dimensions.
existing_points : ndarray or list of lists, default=None
Set of existing points in [0, 1]^d.
random_state : int or RandomState, default=None
Pseudo random number generator state used for random uniform sampling from lists
of possible values instead of scipy.stats distributions.
restarts : int, default=20
Number of random restarts to use for computing the global optimum of the energy
landscape.
Returns
-------
ndarray, shape (n, d)
Set of n points in [0, 1]^d.
Raises
------
ValueError
If the number of existing points exceeds the number of total requested points.
References
----------
.. [1] Steinerberger, S. Dynamically defined sequences with small discrepancy.
Monatsh Math 191, 639–655 (2020). https://doi.org/10.1007/s00605-019-01360-z
"""
random_state = check_random_state(random_state)
if existing_points is None:
X = [random_state.uniform(size=d)]
else:
X = list(existing_points)
if len(X) >= n:
raise ValueError("No more points left to generate.")
n -= len(X)
for _ in range(n):
random_starts = random_state.uniform(size=(restarts, d))
best_value = np.inf
best_point = random_starts[0]
for start in random_starts:
with np.errstate(invalid="ignore"):
result = minimize(
_sb_functional,
x0=start,
bounds=[(0.0, 1.0)] * d,
args=(np.array(X),),
)
if result.fun < best_value:
best_point = result.x
best_value = result.fun
X.append(best_point)
return np.array(X)
def phi(d, n_iter=10):
if d == 1:
return 1.61803398874989484820458683436563
if d == 2:
return 1.32471795724474602596090885447809
x = 2.0000
for _ in range(n_iter):
x = pow(1 + x, 1 / (d + 1))
return x
[docs]def r2_sequence(n, d, seed=0.5):
"""Output ``n`` points of the infinite R2 quasi-random sequence.
Parameters
----------
n : int
Number of points to generate
d : int
Number of dimensions for each point
seed : float in [0, 1], default=0.5
Seed value for the sequence
Returns
-------
z : ndarray, shape (n, d)
``n`` points of the R2 sequence
"""
g = phi(d)
alpha = np.zeros(d)
for j in range(d):
alpha[j] = pow(1 / g, j + 1) % 1
z = np.zeros((n, d))
for i in range(n):
z[i] = (seed + alpha * (i + 1)) % 1
return z