Page MenuHomec4science
Files F60770843

sobol.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 12:25

sobol.py
View Options

# Copyright (C) 2018-2020 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
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['sobolGenerate']
def low2bit(n):
if n < 0:
raise RROMPyException("Only positive integers allowed.")
j = 2
while (n % j) // (j / 2) == 1: j = j * 2
return int(np.log2(j)) - 1
def sobolGenerate(dim_num, n, seed = 0, return_seed = False):
SEEDBITS = 30
v = np.zeros((dim_num, SEEDBITS), dtype = int)
v[:, 0] = 1
v[2:, 1] = np.array([1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3,
1,3,1,1,3,1,3,1,3,1,3])[:max(0,dim_num-2)]
v[3:, 2] = np.array([7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7,
5,1,1,5,7,7,5,1,3,3])[:max(0,dim_num-3)]
v[5:, 3] = np.array([1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1,
11,7,5,15,1,15,11,5,3,1,7,9])[:max(0,dim_num-5)]
v[7:, 4] = np.array([9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5,
23,27,19,21,5,1,17,13,7,15,9,31,9])[:max(0,dim_num-7)]
v[13:, 5] = np.array([37,33,7,5,11,39,63,27,17,15,23,29,3,21,13,31,25,9,49,
33,19,29,11,19,27,15,25])[:max(0,dim_num-13)]
v[19:, 6] = np.array([13,33,115,41,79,17,29,119,75,73,105,7,59,65,21,3,113,
61,89,45,107])[:max(0,dim_num-19)]
v[37:, 7] = np.array([7,23,39])[:max(0,dim_num-37)]
poly = [1,3,7,11,13,19,25,37,59,47,61,55,41,67,97,91,109,103,115,131,193,
137,145,143,241,157,185,167,229,171,213,191,253,203,211,239,247,
285,369,299]
v[0, :] = 1
for i in range(1, dim_num):
j = poly[i]
m = int(1 + np.floor(np.log2(j)))
includ = np.array([int(x) for x in "{0:b}".format(j)], dtype = int)
for j in range(m - 1, SEEDBITS):
newv = v[i, j - m + 1]
l = 1
for k in range(1, m):
l *= 2
if includ[k]:
newv = np.bitwise_xor(newv, l * v[i, j - k])
v[i, j] = newv
v = np.multiply(v, np.power(2, np.arange(SEEDBITS)[::-1]))
recipd = np.power(2., - SEEDBITS)
lastq = np.zeros(dim_num, dtype = int)
####
if seed < 0:
seed = 0
for seed_temp in range(seed):
l = low2bit(seed_temp)
lastq = np.bitwise_xor(lastq, v[:, l])
r = np.empty((n, dim_num))
for seed_temp in range(seed, seed + n):
l = low2bit(seed_temp)
r[seed_temp - seed, :] = lastq * recipd
lastq = np.bitwise_xor(lastq, v[:, l])
if return_seed: return r, seed + n
return r

Event Timeline

c4science · Help