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


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