Fractional factorial sensitivity analysis

SALib/analyze/ff.py

@@ -0,0 +1,84 @@
+'''
+Created on 30 Jun 2015
+
+@author: will2
+'''
+
+from __future__ import print_function
+import numpy as np
+from . import common_args
+from SALib.util import read_param_file, unscale_samples
+from SALib.sample.ff import generate_contrast, sample, extend_bounds
+
+def analyze(problem, X, Y, second_order=False, print_to_console=False):
+
+    problem = extend_bounds(problem)
+    num_vars = problem['num_vars']
+
+    X = generate_contrast(problem)
+
+    main_effect = (1. / (2 * num_vars)) * np.dot(Y, X)
+
+    Si = dict((k, [None] * num_vars)
+              for k in ['names', 'ME'])
+    Si['ME'] = main_effect
+    Si['names'] = problem['names']
+
+    if print_to_console:
+        print("Parameter ME")
+        for j in range(num_vars):
+            print("%s %f" % (problem['names'][j], Si['ME'][j]))
+
+    if second_order == True:
+        interaction_names, interaction_effects = interactions(problem, 
+                                                              Y, 
+                                                              print_to_console)
+
+        Si['names'].append(interaction_names)
+        Si['IE'] = interaction_effects
+
+    return Si
+
+
+def interactions(problem, Y, print_to_console=False):
+    '''
+    Computes the second order effects (interactions) between
+    all combinations of pairs of input factors
+    '''
+
+    names = problem['names']
+    num_vars = problem['num_vars']
+
+    X = generate_contrast(problem)
+
+    ie_names = []
+    IE = []
+
+    for col in range(X.shape[1]):
+        for col_2 in range(col):
+            x = X[:, col] * X[:, col_2]
+            var_names = names[col_2] + names[col]
+            ie_names.append(var_names)
+            IE.append((1. / (2 * num_vars)) * np.dot(Y, x))
+    if print_to_console:
+        [print('%s %f' % (n, i) ) for (n, i) in zip(ie_names, IE) ]
+
+    return ie_names, IE
+
+if __name__ == "__main__":
+
+    parser = common_args.create()
+    parser.add_argument('-X', '--model-input-file', type=str,
+                        required=True, default=None, help='Model input file')
+    parser.add_argument('--max-order', type=int, required=False, default=2,
+                    choices=[1, 2], help='Maximum order of sensitivity indices to calculate')
+    args = parser.parse_args()
+
+    problem = read_param_file(args.paramfile)
+
+    Y = np.loadtxt(args.model_output_file, delimiter=args.delimiter, usecols=(args.column,))
+    X = np.loadtxt(args.model_input_file, delimiter=args.delimiter, ndmin=2)
+    if len(X.shape) == 1:
+        X = X.reshape((len(X), 1))
+    analyze(problem, X, Y, (args.max_order == 2), print_to_console=True)
+

SALib/sample/common_args.py

@@ -5,8 +5,6 @@ def create():
     parser = argparse.ArgumentParser(
         description='Create parameter samples for sensitivity analysis')
 
-    parser.add_argument(
-        '-n', '--samples', type=int, required=True, help='Number of Samples')
     parser.add_argument(
         '-p', '--paramfile', type=str, required=True, help='Parameter Range File')
     parser.add_argument(

SALib/sample/fast_sampler.py

@@ -49,6 +49,9 @@ def sample(problem, N, M=4):
 if __name__ == "__main__":
 
     parser = common_args.create()
+    parser.add_argument(
+        '-n', '--samples', type=int, required=True, help='Number of Samples')
+
     parser.add_argument(
         '-M', type=int, required=False, default=4, help='M coefficient, default 4')
     args = parser.parse_args()

SALib/sample/ff.py

@@ -0,0 +1,70 @@
+'''
+Created on 30 Jun 2015
+
+@author: will2
+'''
+
+from scipy.linalg import hadamard
+import numpy as np
+from . import common_args
+from ..util import scale_samples, read_param_file
+
+
+def find_smallest(num_vars):
+    '''
+    Find the smallest exponent of two that is greater than the number
+    of variables
+    '''
+    for x in range(10):
+        if num_vars <= 2 ** x:
+            return x
+
+
+def extend_bounds(problem):
+
+    num_vars = problem['num_vars']
+    num_ff_vars = 2 ** find_smallest(num_vars)
+    num_dummy_variables = num_ff_vars - num_vars
+
+    bounds = list(problem['bounds'])
+    names = problem['names']
+    if num_dummy_variables > 0:
+        bounds.extend([[0, 1] for x in range(num_dummy_variables)])
+        names.extend(["dummy_" + str(var) for var in range(num_dummy_variables)])
+        problem['bounds'] = bounds
+        problem['names'] = names
+        problem['num_vars'] = num_ff_vars
+
+    return problem
+
+def generate_contrast(problem):
+
+    num_vars = problem['num_vars']
+
+    # Find the smallest n, such that num_vars < k
+    k = [2 ** n for n in range(16)]
+    k_chosen = 2 ** find_smallest(num_vars)
+
+    # Generate the fractional factorial contrast
+    contrast = np.vstack([hadamard(k_chosen), -hadamard(k_chosen)])
+
+    return contrast
+
+def sample(problem):
+
+    contrast = generate_contrast(problem)
+    sample = np.array((contrast + 1.) / 2, dtype=np.float)
+    problem = extend_bounds(problem)
+    scale_samples(sample, problem['bounds'])
+    return sample
+
+if __name__ == "__main__":
+
+    parser = common_args.create()
+    args = parser.parse_args()
+
+    problem = read_param_file(args.paramfile)
+    param_values = sample(problem)
+
+    np.savetxt(args.output, param_values, delimiter=args.delimiter,
+               fmt='%.' + str(args.precision) + 'e')

SALib/sample/finite_diff.py

@@ -47,6 +47,9 @@ def sample(problem, N, delta=0.01):
     parser = common_args.create()
     parser.add_argument('-d', '--delta', type=float, required=False,
                         default=0.01, help='Finite difference step size (percent)')
+    parser.add_argument(
+        '-n', '--samples', type=int, required=True, help='Number of Samples')
+
     args = parser.parse_args()
 
     np.random.seed(args.seed)

SALib/sample/latin.py

@@ -33,7 +33,9 @@ def sample(problem, N):
 
     parser = common_args.create()
     args = parser.parse_args()
-
+    parser.add_argument(
+        '-n', '--samples', type=int, required=True, help='Number of Samples')
+
     np.random.seed(args.seed)
     problem = read_param_file(args.paramfile)
 

SALib/sample/morris.py

@@ -177,6 +177,9 @@ def compute_optimised_trajectories(problem, input_sample, N, k_choices):
 if __name__ == "__main__":
 
     parser = common_args.create()
+
+    parser.add_argument(
+        '-n', '--samples', type=int, required=True, help='Number of Samples')
     parser.add_argument('-l', '--levels', type=int, required=False,
                         default=4, help='Number of grid levels (Morris only)')
     parser.add_argument('--grid-jump', type=int, required=False,

SALib/sample/saltelli.py

@@ -68,6 +68,10 @@ def sample(problem, N, calc_second_order=True):
 if __name__ == "__main__":
 
     parser = common_args.create()
+
+    parser.add_argument(
+        '-n', '--samples', type=int, required=True, help='Number of Samples')
+
     parser.add_argument('--max-order', type=int, required=False, default=2,
                         choices=[1, 2], help='Maximum order of sensitivity indices to calculate')
     args = parser.parse_args()

SALib/test_functions/Sobol_G.py

@@ -14,9 +14,9 @@ def evaluate(values, a=[0, 1, 4.5, 9, 99, 99, 99, 99]):
 
     if type(values) != np.ndarray:
         raise TypeError("The argument `values` must be a numpy ndarray")
-
+
     ltz = np.array(values) < 0
-    gto = np.array(values) > 1
+    gtz = np.array(values) > 1
 
     if ltz.any() == True:
         raise ValueError("Sobol G function called with values less than one")

SALib/tests/test_ff.py

@@ -0,0 +1,170 @@
+'''
+Created on 30 Jun 2015
+
+@author: will2
+'''
+import numpy as np
+from numpy.testing import assert_equal, assert_allclose
+from ..sample.ff import sample, find_smallest, extend_bounds
+from ..analyze.ff import analyze, interactions
+
+
+def test_find_smallest():
+    '''
+    '''
+    num_vars = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16, 17, 31, 32, 33]
+    expected = [0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6]
+    for x, y in zip(num_vars, expected):
+        actual = find_smallest(x)
+        assert_equal(actual, y)
+
+
+def test_extend_bounds():
+    problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T,
+           'num_vars': 12,
+           'names': ["x" + str(x + 1) for x in range(12)]
+           }
+    actual = extend_bounds(problem)
+    expected = {'names': ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'dummy_0', 'dummy_1', 'dummy_2', 'dummy_3'], 
+                'bounds': [np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([-1,  1]), np.array([-1,  1]), 
+                           np.array([0, 1]), np.array([0, 1]), 
+                           np.array([0, 1]), np.array([0, 1])], 
+                'num_vars': 16}
+
+    assert_equal(actual, expected)
+
+def test_ff_sample():
+    problem = {'bounds': [[0., 1.], [0., 1.], [0., 1.], [0., 1.]],
+               'num_vars': 4,
+               'names': ['x1', 'x2', 'x3', 'x4']}
+    actual = sample(problem)
+    expected = np.array([[ 1, 1, 1, 1],
+                         [ 1, 0, 1, 0],
+                         [ 1, 1, 0, 0],
+                         [ 1, 0, 0, 1],
+                         [0, 0, 0, 0],
+                         [0, 1, 0, 1],
+                         [0, 0, 1, 1],
+                         [0, 1, 1, 0]], dtype=np.float)
+    assert_equal(actual, expected)
+
+
+def test_ff_sample_scaled():
+    '''
+    '''
+    problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]],
+               'num_vars': 4,
+               'names': ['x1', 'x2', 'x3', 'x4']}
+    actual = sample(problem)
+    expected = np.array([[ 2.5, 1, 1, 1],
+                         [ 2.5, 0, 1, 0],
+                         [ 2.5, 1, 0, 0],
+                         [ 2.5, 0, 0, 1],
+                         [0, 0, 0, 0],
+                         [0, 1, 0, 1],
+                         [0, 0, 1, 1],
+                         [0, 1, 1, 0]], dtype=np.float)
+    assert_equal(actual, expected)
+
+
+def test_ff_analyze():
+    '''
+    '''
+
+    problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]],
+               'num_vars': 4,
+               'names': ['x1', 'x2', 'x3', 'x4']}
+    X = np.array([[ 1, 1, 1, 1],
+                  [ 1, 0, 1, 0],
+                  [ 1, 1, 0, 0],
+                  [ 1, 0, 0, 1],
+                  [0, 0, 0, 0],
+                  [0, 1, 0, 1],
+                  [0, 0, 1, 1],
+                  [0, 1, 1, 0]], dtype=np.float)
+    Y = np.array([1.5, 1, 1.5, 1, 2, 2.5, 2, 2.5], dtype=np.float)
+    actual = analyze(problem, X, Y)
+    expected = {'ME': np.array([ -0.5 ,  0.25,  0.  ,  0.  ]), 'names': ['x1', 'x2', 'x3', 'x4']}
+    assert_equal(actual, expected)
+
+
+def test_ff_example():
+    '''
+    '''
+
+    problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T,
+               'num_vars': 12,
+               'names': ["x" + str(x + 1) for x in range(12)]
+               }
+
+    X = sample(problem)
+    Y = X[:, 0] + 2 * X[:, 1] + 3 * X[:, 2] + 4 * X[:, 6] * X[:, 11]
+
+    expected = np.array([10, -2, 4, -8, 2, 6, -4,
+                         0, 2, 6, -4, 0, 10, -2, 4, -8,
+                         - 2, -6, 4, 0, -10, 2, -4, 8,
+                          - 10, 2, -4, 8, -2, -6, 4, 0])
+
+    assert_equal(Y, expected)
+
+    Si = analyze(problem, X, Y)
+
+    expected = np.array([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=np.float)
+    assert_equal(expected, Si['ME'])
+
+
+def test_interactions_from_saltelli():
+    '''
+    '''
+    problem = {'bounds': np.repeat([-1, 1], 12).reshape(2, 12).T,
+               'num_vars': 12,
+               'names': ["x" + str(x + 1) for x in range(12)]
+               }
+
+    X = sample(problem)
+
+    Y = np.array([10, -2, 4, -8, 2, 6, -4, 0,
+                   2, 6, -4, 0, 10, -2, 4, -8,
+                  - 2, -6, 4, 0, -10, 2, -4, 8,
+                 - 10, 2, -4, 8, -2, -6, 4, 0])
+
+    Si = analyze(problem, X, Y, second_order=True)
+    actual = Si['IE']
+    expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+    assert_equal(actual, expected)
+
+
+def test_interactions():
+    '''
+    '''
+    problem = {'bounds': [[0., 2.5], [0., 1.], [0., 1.], [0., 1.]],
+               'num_vars': 4,
+               'names': ['x1', 'x2', 'x3', 'x4']}
+    X = np.array([[ 2.5, 1.0, 1.0, 1.0],
+                  [ 2.5, 0, 1.0, 0],
+                  [ 2.5, 1.0, 0, 0],
+                  [ 2.5, 0, 0, 1.0],
+                  [0, 0, 0, 0],
+                  [0, 1.0, 0, 1.0],
+                  [0, 0, 1.0, 1.0],
+                  [0, 1.0, 1.0, 0]], dtype=np.float)
+    Y = X[:, 0] + (0.1 * X[:, 1]) + ((1.2 * X[:, 2]) * (1.3 + X[:, 3]))
+#     Y = np.array([1.5, 1, 1.5, 1, 2, 2.5, 2, 2.5], dtype=np.float)
+    ie_names, ie = interactions(problem, Y, print_to_console=True)
+    actual = ie
+    assert_allclose(actual, [0.3, 0, 0, 0, 0, 0.3], rtol=1e-4, atol=1e-4)