Merge pull request SALib#138 from SALib/bugfix-136-sobol-sequence

fix for sobol sequence: first row should be zero not empty. add test …

SALib/analyze/delta.py

@@ -96,7 +96,7 @@ def calc_delta(Y, Ygrid, X, m):
 # Plischke et al. 2013 bias reduction technique (eqn 30)
 
 def bias_reduced_delta(Y, Ygrid, X, m, num_resamples, conf_level):
-    d = np.empty(num_resamples)
+    d = np.zeros(num_resamples)
     d_hat = calc_delta(Y, Ygrid, X, m)
 
     for i in range(num_resamples):
@@ -119,7 +119,7 @@ def sobol_first(Y, X, m):
 
 
 def sobol_first_conf(Y, X, m, num_resamples, conf_level):
-    s = np.empty(num_resamples)
+    s = np.zeros(num_resamples)
 
     for i in range(num_resamples):
         r = np.random.randint(len(Y), size=len(Y))

SALib/analyze/dgsm.py

@@ -51,10 +51,10 @@ def analyze(problem, X, Y, num_resamples=1000,
     if not 0 < conf_level < 1:
         raise RuntimeError("Confidence level must be between 0-1.")
 
-    base = np.empty(N)
-    X_base = np.empty((N, D))
-    perturbed = np.empty((N, D))
-    X_perturbed = np.empty((N, D))
+    base = np.zeros(N)
+    X_base = np.zeros((N, D))
+    perturbed = np.zeros((N, D))
+    X_perturbed = np.zeros((N, D))
     step = D + 1
 
     base = Y[0:Y.size:step]
@@ -65,7 +65,7 @@ def analyze(problem, X, Y, num_resamples=1000,
 
     # First order (+conf.) and Total order (+conf.)
     keys = ('vi', 'vi_std', 'dgsm', 'dgsm_conf')
-    S = dict((k, np.empty(D)) for k in keys)
+    S = dict((k, np.zeros(D)) for k in keys)
     if print_to_console:
         print("Parameter %s %s %s %s" % keys)
 
@@ -98,7 +98,7 @@ def calc_dgsm(base, perturbed, x_delta, bounds, num_resamples, conf_level):
     vi, _ = calc_vi(base, perturbed, x_delta)
     dgsm = vi * (bounds[1] - bounds[0]) ** 2 / (D * np.pi ** 2)
 
-    s = np.empty(num_resamples)
+    s = np.zeros(num_resamples)
     for i in range(num_resamples):
         r = np.random.randint(len(base), size=len(base))
         s[i], _ = calc_vi(base[r], perturbed[r], x_delta[r])

SALib/analyze/fast.py

@@ -60,7 +60,7 @@ def analyze(problem, Y, M=4, print_to_console=False):
         exit()
 
     # Recreate the vector omega used in the sampling
-    omega = np.empty([D])
+    omega = np.zeros([D])
     omega[0] = math.floor((N - 1) / (2 * M))
     m = math.floor(omega[0] / (2 * M))
 

SALib/analyze/morris.py

@@ -159,7 +159,7 @@ def compute_grouped_sigma(ungrouped_sigma, group_matrix):
     sigma_masked = np.ma.masked_array(ungrouped_sigma * group_matrix.T,
                                         mask=(group_matrix^1).T)
     sigma_agg = np.ma.mean(sigma_masked, axis=1)
-    sigma = np.empty(group_matrix.shape[1], dtype=np.float)
+    sigma = np.zeros(group_matrix.shape[1], dtype=np.float)
     np.copyto(sigma, sigma_agg, where=group_matrix.sum(axis=0)==1 )
     np.copyto(sigma, np.NAN, where=group_matrix.sum(axis=0)!=1 )
 
@@ -241,8 +241,8 @@ def compute_mu_star_confidence(ee, num_trajectories, num_resamples, conf_level):
     to produce a histogram of resampled mu_star metrics.
     This resample is used to produce a confidence interval.
     '''
-    ee_resampled = np.empty([num_trajectories])
-    mu_star_resampled = np.empty([num_resamples])
+    ee_resampled = np.zeros([num_trajectories])
+    mu_star_resampled = np.zeros([num_resamples])
 
     if not 0 < conf_level < 1:
         raise ValueError("Confidence level must be between 0-1.")

SALib/analyze/sobol.py

@@ -149,20 +149,20 @@ def second_order(A, ABj, ABk, BAj, B):
 
 def create_Si_dict(D, calc_second_order):
     # initialize empty dict to store sensitivity indices
-    S = dict((k, np.empty(D)) for k in ('S1', 'S1_conf', 'ST', 'ST_conf'))
+    S = dict((k, np.zeros(D)) for k in ('S1', 'S1_conf', 'ST', 'ST_conf'))
 
     if calc_second_order:
-        S['S2'] = np.empty((D, D))
+        S['S2'] = np.zeros((D, D))
         S['S2'][:] = np.nan
-        S['S2_conf'] = np.empty((D, D))
+        S['S2_conf'] = np.zeros((D, D))
         S['S2_conf'][:] = np.nan
 
     return S
 
 
 def separate_output_values(Y, D, N, calc_second_order):
-    AB = np.empty((N, D))
-    BA = np.empty((N, D)) if calc_second_order else None
+    AB = np.zeros((N, D))
+    BA = np.zeros((N, D)) if calc_second_order else None
     step = 2 * D + 2 if calc_second_order else D + 2
 
     A = Y[0:Y.size:step]

SALib/sample/fast_sampler.py

@@ -33,7 +33,7 @@ def sample(problem, N, M=4):
 
     D = problem['num_vars']
 
-    omega = np.empty([D])
+    omega = np.zeros([D])
     omega[0] = math.floor((N - 1) / (2 * M))
     m = math.floor(omega[0] / (2 * M))
 
@@ -46,8 +46,8 @@ def sample(problem, N, M=4):
     s = (2 * math.pi / N) * np.arange(N)
 
     # Transformation to get points in the X space
-    X = np.empty([N * D, D])
-    omega2 = np.empty([D])
+    X = np.zeros([N * D, D])
+    omega2 = np.zeros([D])
 
     for i in range(D):
         omega2[i] = omega[0]

SALib/sample/finite_diff.py

@@ -20,7 +20,7 @@ def sample(problem, N, delta=0.01):
     base_sequence = sobol_sequence.sample(N + skip_values, D)
     # scale before finite differencing
     scale_samples(base_sequence, problem['bounds'])
-    dgsm_sequence = np.empty([N * (D + 1), D])
+    dgsm_sequence = np.zeros([N * (D + 1), D])
 
     index = 0
 

SALib/sample/latin.py

@@ -24,8 +24,8 @@ def sample(problem, N):
     """
     D = problem['num_vars']
 
-    result = np.empty([N, D])
-    temp = np.empty([N])
+    result = np.zeros([N, D])
+    temp = np.zeros([N])
     d = 1.0 / N
 
     for i in range(D):

SALib/sample/morris.py

@@ -107,7 +107,7 @@ def sample_oat(problem, N, num_levels, grid_jump):
 
     # grid step delta, and final sample matrix X
     delta = grid_jump / (num_levels - 1)
-    X = np.empty([N * (D + 1), D])
+    X = np.zeros([N * (D + 1), D])
 
     # Create N trajectories. Each trajectory contains D+1 parameter sets.
     # (Starts at a base point, and then changes one parameter at a time)
@@ -123,7 +123,7 @@ def sample_oat(problem, N, num_levels, grid_jump):
             P[i, perm[i]] = 1
 
         # starting point for this trajectory
-        x_base = np.empty([D + 1, D])
+        x_base = np.zeros([D + 1, D])
         for i in range(D):
             x_base[:, i] = (
                 rd.choice(np.arange(num_levels - grid_jump))) / (num_levels - 1)
@@ -152,7 +152,7 @@ def sample_groups(problem, N, num_levels, grid_jump):
 
     k = G.shape[0]
     g = G.shape[1]
-    sample = np.empty((N * (g + 1), k))
+    sample = np.zeros((N * (g + 1), k))
     sample = np.array([generate_trajectory(G, num_levels, grid_jump) for n in range(N)])
     return sample.reshape((N * (g + 1), k))
 

SALib/sample/morris_util.py

@@ -72,7 +72,7 @@ def generate_x_star(k, num_levels, grid_step):
     Generate an 1-by-k array to represent initial position for EE
     This should be a randomly generated array in the p level grid :math:\omega
     '''
-    x_star = np.empty(k)
+    x_star = np.zeros(k)
 
     delta = compute_delta(num_levels)
     bound = 1 - delta
@@ -171,7 +171,7 @@ def find_most_distant(input_sample, N, num_params, k_choices, groups=None):
     # Generate a list of all the possible combinations
     # combos = np.array([x for x in combinations(range(N),k_choices)])
     combo_gen = combinations(list(range(N)), k_choices)
-    scores = np.empty(number_of_combinations, dtype=np.float32)
+    scores = np.zeros(number_of_combinations, dtype=np.float32)
     # Generate the pairwise indices once
     pairwise = np.array([y for y in combinations(list(range(k_choices)), 2)])
 

SALib/sample/saltelli.py

@@ -43,9 +43,9 @@ def sample(problem, N, calc_second_order=True):
     base_sequence = sobol_sequence.sample(N + skip_values, 2 * D)
 
     if calc_second_order:
-        saltelli_sequence = np.empty([(2 * Dg + 2) * N, D])
+        saltelli_sequence = np.zeros([(2 * Dg + 2) * N, D])
     else:
-        saltelli_sequence = np.empty([(Dg + 2) * N, D])
+        saltelli_sequence = np.zeros([(Dg + 2) * N, D])
     index = 0
 
     for i in range(skip_values, N + skip_values):

SALib/sample/sobol_sequence.py

@@ -51,8 +51,8 @@
 def sample(N, D):
     """Generate (N x D) numpy array of Sobol sequence samples"""
     scale = 31
-    result = np.empty([N, D])
-
+    result = np.zeros([N, D])
+    
     if D > len(directions) + 1:
         raise ValueError("Error in Sobol sequence: not enough dimensions")
 
@@ -62,7 +62,7 @@ def sample(N, D):
         raise ValueError("Error in Sobol sequence: not enough bits")
 
     for i in range(D):
-        V = np.empty(L + 1, dtype=long)
+        V = np.zeros(L + 1, dtype=long)
 
         if i == 0:
             for j in range(1, L + 1):

SALib/test_functions/Ishigami.py

@@ -11,7 +11,7 @@
 # x2: 0.4424
 # x3: 0.0
 def evaluate(values):
-    Y = np.empty([values.shape[0]])
+    Y = np.zeros([values.shape[0]])
     A = 7
     B = 0.1
 

SALib/test_functions/Sobol_G.py

@@ -25,7 +25,7 @@ def evaluate(values, a=None):
     elif gto.any() == True:
         raise ValueError("Sobol G function called with values greater than one")
 
-    Y = np.empty([values.shape[0]])
+    Y = np.zeros([values.shape[0]])
 
     for i, row in enumerate(values):
         Y[i] = 1.0

docs/basics.rst

@@ -121,7 +121,7 @@ will loop over each sample input and evaluate the model:
 
 .. code:: python
 
-    Y = np.empty([param_values.shape[0]])
+    Y = np.zeros([param_values.shape[0]])
 
     for i, X in enumerate(param_values):
         Y[i] = evaluate_model(X)

tests/test_morris_util.py

@@ -253,7 +253,7 @@ def test_find_maximum():
 
 
 def test_catch_combos_too_large():
-    N = 1e6
+    N = int(1e6)
     k_choices = 4
     num_params = 2
     input_sample = np.random.random_sample((N, num_params))

tests/test_sobol.py

@@ -6,7 +6,7 @@
 from scipy.stats import norm
 
 from SALib.analyze import sobol
-from SALib.sample import saltelli
+from SALib.sample import saltelli, sobol_sequence
 from SALib.test_functions import Ishigami, Sobol_G
 from SALib.util import read_param_file
 
@@ -18,6 +18,15 @@ def setup_samples(N = 500, calc_second_order = True):
     return problem,param_values
 
 
+def test_sobol_sequence():
+    # example from Joe & Kuo: http://web.maths.unsw.edu.au/~fkuo/sobol/
+    S = sobol_sequence.sample(10,3)
+    expected = [[0,0,0],[0.5,0.5,0.5],[0.75,0.25,0.25],[0.25,0.75,0.75],
+                [0.375,0.375,0.625],[0.875,0.875,0.125],[0.625,0.125,0.875],
+                [0.125,0.625,0.375],[0.1875,0.3125,0.9375],[0.6875,0.8125,0.4375]]
+    assert_allclose(S, expected, atol=5e-2, rtol=1e-1)
+
+
 def test_sample_size_second_order():
     N = 500
     D = 3