Page MenuHomec4science
Files F60115181

halton.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 15:04

halton.py
View Options

# 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']
haltonPrimes = [2, 3]
def haltonAddPrimes(n):
global haltonPrimes
for j in range(len(haltonPrimes), n):
primeBase = haltonPrimes[- 1] + 2
while True:
for i in range(1, j):
if primeBase % haltonPrimes[i] == 0: break
else:
break
primeBase += 2
haltonPrimes += [primeBase]
return haltonPrimes
def haltonGenerate(dim, size, seed = 0, step = True, return_seed = False):
haltonPrimes = haltonAddPrimes(dim * (1 + 1 * step) + 1)
base = np.array(haltonPrimes[: dim])
stepSize = haltonPrimes[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
if return_seed: return seq, stepSize * size + seed
return seq

Event Timeline

c4science · Help