# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#

import numpy as np

__all__ = ['haltonGenerate']

def haltonGenerate(dim, size, seed = 0, step = True):
    seq = np.empty((size, dim), dtype = float)
    primes = [2, 3]
    for j in range(2, dim * (1 + 1 * step) + 1):
        primeBase = primes[j - 1] + 2
        while True:
            for i in range(1, j):
                if primeBase % primes[i] == 0: break
            else:
                break
            primeBase += 2
        primes += [primeBase]
    base = np.array(primes[: dim])
    stepSize = primes[dim * 2 - 1] if step else 1
    seq = np.zeros((size, dim), dtype = float)
    seqBase = np.tile((np.arange(seed, stepSize * size + seed, stepSize)
                     + 1).reshape(-1, 1), [1, dim])
    powerB = 1
    while seqBase[-1, 0] > 0:
        seqBase, remainder = np.divmod(seqBase, base)
        seq += remainder / np.power(base, powerB)
        powerB += 1
    return seq