Repository: salib/salib
Issue: Morris: Identical trajectories cause differences between optimal trajectories

If identical trajectories were generated (which is only likely with very small numbers of parameters), then duplicate entries can be made in the distance matrix which leads to duplicate (but identical) combinations of trajectories than can be chosen.


Changed files:
SALib/sample/morris_util.py, SALib/sample/optimal_trajectories.py, SALib/tests/test_optimal_trajectories.py, SALib/tests/test_morris_util.py, SALib/analyze/morris.py

diff --git a/SALib/analyze/morris.py b/SALib/analyze/morris.py
index c83eace..ed3b566 100755
--- a/SALib/analyze/morris.py
+++ b/SALib/analyze/morris.py
@@ -61,7 +61,7 @@ def analyze(problem, X, Y,
                     'sigma'][j], Si['mu_star'][j], Si['mu_star_conf'][j]))
         return Si
     elif groups is not None:
-        Si['sigma'] = np.std(np.abs(ee), axis=1, ddof=1)
+        Si['sigma'] = np.std(ee, axis=1, ddof=1)
         # if there are groups, then the elementary effects returned need to be
         # computed over the groups of variables, rather than the individual variables
         Si_grouped = dict((k, [None] * num_vars)
@@ -86,23 +86,35 @@ def analyze(problem, X, Y,
 
 
 def compute_grouped_metric(ungrouped_metric, group_matrix):
+    '''
+    Computes the mean value for the groups of parameter values in the
+    argument ungrouped_metric
+    '''
+
+    group_matrix = np.array(group_matrix, dtype=np.bool)
 
-    group_matrix = np.array(group_matrix)
-    grouped_metric = np.divide(np.dot(ungrouped_metric, group_matrix), np.sum(group_matrix, 0))
+    mu_star_masked = np.ma.masked_array(ungrouped_metric * group_matrix.T, 
+                                        mask=(group_matrix^1).T)
+    mean_of_mu_star = np.ma.mean(mu_star_masked, axis=1)
 
-    return grouped_metric.T
+    return mean_of_mu_star
 
 
 def compute_grouped_sigma(Si, group_matrix):
+    '''
+    Computes a pooled standard deviation from the groups of parameters defined
+    but the group_matrix, from the sum of the mean of variances and variance
+    of the means of the (absolute) elementary effects.
+    '''
     
     group_matrix = np.array(group_matrix, dtype=np.bool)
 
-    mu_star_masked = np.ma.masked_array(Si['mu_star'] * group_matrix.T, mask=(group_matrix^1).T)
-    sigma_masked = np.ma.masked_array(Si['sigma'] * group_matrix.T, mask=(group_matrix^1).T)**2
+    mu_star_masked = np.ma.masked_array(Si['mu_star'] * group_matrix.T, 
+                                        mask=(group_matrix^1).T)
+    sigma_masked = np.ma.masked_array(Si['sigma'] * group_matrix.T, 
+                                      mask=(group_matrix^1).T)**2
     
     mean_of_variances_of_ee = np.ma.mean(sigma_masked, axis=1)
-#     variance_of_means = np.var(Si['mu_star'] * group_matrix.T, axis=1)
-# ``var = mean(abs(x - x.mean())**2)``
 
     variance_of_means_of_ee = np.ma.var(mu_star_masked, axis=1, ddof=1)
 
diff --git a/SALib/sample/morris_util.py b/SALib/sample/morris_util.py
index 489d7d5..944edbf 100644
--- a/SALib/sample/morris_util.py
+++ b/SALib/sample/morris_util.py
@@ -100,7 +100,7 @@ def check_input_sample(input_sample, num_params, N):
     "Input sample does not match number of parameters or groups"
     assert np.any((input_sample >= 0) | (input_sample <= 1)), \\
     "Input sample must be scaled between 0 and 1"
-
+            
 
 def compute_distance(m, l):
     '''
@@ -109,8 +109,11 @@ def compute_distance(m, l):
 
     if np.shape(m) != np.shape(l):
         raise ValueError("Input matrices are different sizes")
-
-    distance = np.array(np.sum(cdist(m, l)), dtype=np.float32)
+    if np.array_equal(m, l):
+        print("Trajectory %s and %s are equal" % (m, l))
+        distance = 0
+    else:
+        distance = np.array(np.sum(cdist(m, l)), dtype=np.float32)
 
     return distance
 
@@ -125,6 +128,7 @@ def compute_distance_matrix(input_sample, N, num_params, groups=None):
         check_input_sample(input_sample, num_groups, N)
     else:
         check_input_sample(input_sample, num_params, N)
+    
     index_list = make_index_list(N, num_params, num_groups)
     distance_matrix = np.zeros((N, N), dtype=np.float32)
 
@@ -132,7 +136,7 @@ def compute_distance_matrix(input_sample, N, num_params, groups=None):
         input_1 = input_sample[index_list[j]]
         for k in range(j + 1, N):
             input_2 = input_sample[index_list[k]]
-            distance_matrix[k, j] = compute_distance(input_1, input_2)
+            distance_matrix[k, j] = compute_distance(input_1, input_2)    
     return distance_matrix
 
 
@@ -254,7 +258,11 @@ def find_optimum_combination(input_sample, N, num_params, k_choices, groups=None
                                k_choices,
                                groups)
 
+    print(scores)
+
     maximum_combo = find_maximum(scores, N, k_choices)
 
+    print(maximum_combo)
+
     return maximum_combo
 
diff --git a/SALib/sample/optimal_trajectories.py b/SALib/sample/optimal_trajectories.py
index 1ceebb2..a702033 100644
--- a/SALib/sample/optimal_trajectories.py
+++ b/SALib/sample/optimal_trajectories.py
@@ -35,7 +35,7 @@ def model(N, k_choices, distance_matrix):
     m = Model("distance1")
     I = range(N)
 
-    distance_matrix = distance_matrix / distance_matrix.max()
+    distance_matrix = np.array(distance_matrix / distance_matrix.max(), dtype=np.float64)
 
     dm = distance_matrix ** 2
 
@@ -76,7 +76,7 @@ def return_max_combo(input_data, N, num_params, k_choices, groups=None):
     m.params.Threads = 0
     # m.params.MIPGap=0.03
 
-    m.write("model.lp")
+#     m.write("model.lp")
     m.ModelSense = GRB.MAXIMIZE
     m.optimize()
     if m.Status == GRB.OPTIMAL:
@@ -89,6 +89,9 @@ def return_max_combo(input_data, N, num_params, k_choices, groups=None):
                 x_vars.append(v.VarName)
     b = [re.findall("\\d{1,}", str(v)) for v in x_vars]
     maximum_combo = list(set([int(y) for z in b for y in z]))
+    
+    print(maximum_combo)
+    
     return sorted(maximum_combo)
 
 
diff --git a/SALib/tests/test_morris_util.py b/SALib/tests/test_morris_util.py
index fb9eace..385be3b 100644
--- a/SALib/tests/test_morris_util.py
+++ b/SALib/tests/test_morris_util.py
@@ -1,12 +1,13 @@
 from __future__ import division
 
 from nose import with_setup
-from nose.tools import assert_almost_equal, raises
-from numpy.testing import assert_equal
+from nose.tools import raises
+from numpy.testing import assert_equal, assert_array_equal, \\
+                          assert_almost_equal, assert_allclose
 
 import numpy as np
 
-from ..sample.morris_util import generate_P_star, \\
+from SALib.sample.morris_util import generate_P_star, \\
                                 compute_B_star, \\
                                 compute_delta, \\
                                 generate_trajectory, \\
@@ -47,7 +48,7 @@ def test_compute_delta():
     output = [compute_delta(f) for f in fixture]
     desired = np.array([1.00, 0.75, 0.66, 0.62,
                         0.60, 0.58, 0.57, 0.56])
-    np.testing.assert_almost_equal(output, desired, decimal=2)
+    assert_almost_equal(output, desired, decimal=2)
 
 
 def test_generate_trajectory():
@@ -83,7 +84,7 @@ def test_compute_B_star():
     desired = np.array([[1. / 3, 1, 0], [1, 1, 0], [1, 1. / 3, 2. / 3]])
 
     output = compute_B_star(J, x_star, delta, B, G, P_star, D_star)
-    np.testing.assert_array_equal(output, desired)
+    assert_array_equal(output, desired)
 
 
 def test_distance():
@@ -93,7 +94,17 @@ def test_distance():
     input_1 = np.matrix([[0, 1 / 3.], [0, 1.], [2 / 3., 1.]], dtype=np.float32)
     input_3 = np.matrix([[2 / 3., 0], [2 / 3., 2 / 3.], [0, 2 / 3.]], dtype=np.float32)
     output = compute_distance(input_1, input_3)
-    assert_almost_equal(output, 6.18, places=2)
+    assert_allclose(output, 6.18, atol=1e-2)
+
+
+def test_distance_of_identical_matrices_is_min():
+    input_1 = np.matrix([[ 1.          ,1.        ],
+                         [ 1.          ,0.33333333],
+                         [ 0.33333333  ,0.33333333]])
+    input_2 = input_1.copy()
+    actual = compute_distance(input_1, input_2)
+    desired = 0
+    assert_allclose(actual, desired, atol=1e-2)
 
 
 def test_distance_fail_with_difference_size_ip():
@@ -121,7 +132,7 @@ def test_compute_distance_matrix():
     expected[3, :] = [6.89, 6.18, 6.57, 0, 0, 0]
     expected[4, :] = [6.18, 5.31, 5.41, 5.5, 0, 0]
     expected[5, :] = [7.52, 5.99, 5.52, 7.31, 5.77, 0]
-    np.testing.assert_allclose(output, expected, rtol=1e-2)
+    assert_allclose(output, expected, rtol=1e-2)
 
 
 def test_combo_from_find_most_distant():
@@ -137,7 +148,7 @@ def test_combo_from_find_most_distant():
     scores = find_most_distant(sample_inputs, N, num_params, k_choices)
     output = find_maximum(scores, N, k_choices)
     expected = [0, 2, 3, 5]  # trajectories 1, 3, 4, 6
-    np.testing.assert_equal(output, expected)
+    assert_equal(output, expected)
 
 
 def test_scores_from_find_most_distant():
@@ -157,7 +168,7 @@ def test_scores_from_find_most_distant():
                 15.685, 16.098, 14.049, 15.146, 14.333, 14.807, 14.825],
                 dtype=np.float32)
 
-    np.testing.assert_allclose(output, expected, rtol=1e-1, atol=1e-2)
+    assert_allclose(output, expected, rtol=1e-1, atol=1e-2)
 
 
 def test_find_maximum():
diff --git a/SALib/tests/test_optimal_trajectories.py b/SALib/tests/test_optimal_trajectories.py
index 3c96be9..0a51511 100644
--- a/SALib/tests/test_optimal_trajectories.py
+++ b/SALib/tests/test_optimal_trajectories.py
@@ -1,14 +1,15 @@
 from unittest import skipUnless
 
+import numpy as np
 from nose.tools import raises, with_setup
 from numpy.testing import assert_equal
 
-from ..sample.morris import sample_oat, \\
+from SALib.sample.morris import sample_oat, \\
                             find_optimum_combination, \\
                             compute_optimised_trajectories, \\
                             sample_groups
-from ..sample.optimal_trajectories import return_max_combo
-from ..util import read_param_file
+from SALib.sample.optimal_trajectories import return_max_combo
+from SALib.util import read_param_file
 from .test_util import setup_function
 
 
@@ -171,20 +172,41 @@ def test_optimised_trajectories_without_groups():
     """
     Tests that the optimisation problem gives
     the same answer as the brute force problem
-    (for small values of `k_choices` and `N`)
+    (for small values of `k_choices` and `N`),
+    particularly when there are two or more identical
+    trajectories
     """
     
     N = 6
     param_file = "SALib/tests/test_params.txt"
     problem = read_param_file(param_file)
     num_levels = 4
-    grid_jump = num_levels / 2
     k_choices = 4
 
     num_params = problem['num_vars']
     groups = problem['groups']
 
-    input_sample = sample_oat(problem, N, num_levels, grid_jump)
+    # 6 trajectories, with 5th and 6th identical
+    input_sample = np.array([[ 0.33333333,  0.66666667],
+                             [ 1.          ,0.66666667],
+                             [ 1.          ,0.        ],
+                             [ 0.          ,0.33333333],
+                             [ 0.          ,1.        ],
+                             [ 0.66666667  ,1.        ],
+                             [ 0.66666667  ,0.33333333],
+                             [ 0.66666667  ,1.        ],
+                             [ 0.          ,1.        ],
+                             [ 0.66666667  ,1.        ],
+                             [ 0.66666667  ,0.33333333],
+                             [ 0.          ,0.33333333],
+                             [ 1.          ,1.        ],
+                             [ 1.          ,0.33333333],
+                             [ 0.33333333  ,0.33333333],
+                             [ 1.          ,1.        ],
+                             [ 1.          ,0.33333333],
+                             [ 0.33333333  ,0.33333333]], dtype=np.float32)
+
+    print(input_sample)
 
     # From gurobi optimal trajectories     
     actual = return_max_combo(input_sample,
@@ -198,6 +220,7 @@ def test_optimised_trajectories_without_groups():
                                        num_params,
                                        k_choices,
                                        groups)
+    
     assert_equal(actual, desired)