Merge 9febe11 into c1de4e4

README.md

@@ -20,25 +20,22 @@ Python implementations of commonly used sensitivity analysis methods. Useful in
 
 ### Create a parameter file
 
-To get started, create a file describing the sampling ranges for the parameters in the model. Parameter files should be created with 3 columns: `[name, lower bound, upper bound]`:
+To get started, create a file describing the sampling ranges for the parameters in the model. Parameter files should be created with 3 columns:
 ```
+# name lower_bound upper_bound
 P1 0.0 1.0
 P2 0.0 5.0
 ...etc.
 ```
-Lines beginning with `#` will be treated as comments and ignored.
-
-### Or if using Morris... 
-
-You can create a parameter file with groups.
+Lines beginning with `#` will be treated as comments and ignored. The Morris method also supports groups of input factors, which can be specified with a fourth column:
 ```
-#Parameter,Lower Bound, Upper Bound, Group Name
-P 1,0.0,1.0,Group 1
-P 2,0.0,5.0,Group 2
-P 3,0.0,5.0,Group 3
+# name lower_bound upper_bound group_name
+P1 0.0 1.0 Group_1
+P2 0.0 5.0 Group_2
+P3 0.0 5.0 Group_2
 ...etc.
 ```
-You can use delimiters if you wish to include whitespace in your parameter or group names.
+Parameter files can also be comma-delimited if your parameter names or group names contain spaces. This should be detected automatically.
 
 ### Generate samples
 
@@ -50,7 +47,7 @@ python -m SALib.sample.saltelli \
      -o model_input.txt \
 ```
 
-Other methods include `SALib.sample.morris` and `SALib.sample.fast_sampler`. For an explanation of all command line options, [see the examples here](https://github.com/jdherman/SALib/tree/master/examples).
+Other methods include `SALib.sample.morris` and `SALib.sample.fast_sampler`. For an explanation of all command line options for each method, [see the examples here](https://github.com/jdherman/SALib/tree/master/examples).
 
 Or, generate samples from Python:
 ```python
@@ -94,7 +91,7 @@ Si = sobol.analyze(param_file, 'model_output.txt', column = 0, print_to_console=
 # e.g. Si['S1'] contains the first-order index for each parameter, in the same order as the parameter file
 ```
 
-Check out the [examples](https://github.com/jdherman/SALib/tree/master/examples) for a full description of command line and keyword options for all of the methods.
+Check out the [examples](https://github.com/jdherman/SALib/tree/master/examples) for a full description of options for each method.
 
 
 ### License

SALib/analyze/morris.py

@@ -28,7 +28,7 @@ def analyze(pfile,
 
     delta = grid_jump / (num_levels - 1)
 
-    param_file = read_param_file(pfile,True)
+    param_file = read_param_file(pfile)
     Y = np.loadtxt(output_file, delimiter=delim, usecols=(column,))
     X = np.loadtxt(input_file, delimiter=delim, ndmin=2)
     if len(X.shape) == 1:

SALib/sample/morris.py

@@ -3,9 +3,7 @@
 import random as rd
 from . import common_args
 from . sample import Sample
-from ..sample import morris_oat, morris_groups, morris_optimal
-from ..util import read_param_file
-from collections import Iterable
+from . morris_util import *
 
 
 class Morris(Sample):
@@ -36,84 +34,113 @@ def __init__(self, parameter_file, samples, \
                  num_levels, grid_jump, \
                  optimal_trajectories=None):
 
-        self.groups = None
-        self.group_names = None
-        self.parameter_file = parameter_file
-        self.samples = samples
+        Sample.__init__(self, parameter_file, samples)
         self.num_levels = num_levels
-        self.grid_jump = grid_jump
-        pf = read_param_file(self.parameter_file, True)
-        self.num_vars = pf['num_vars']
-        self.bounds = pf['bounds']
-        self.parameter_names = pf['names']
-        if pf['groups'] is not None:
-            self.groups, self.group_names = pf['groups']
+        self.grid_jump = grid_jump        
         self.optimal_trajectories = optimal_trajectories
 
-        if self.optimal_trajectories != None:
-            # Check to ensure that fewer optimal trajectories than samples are
-            # requested, otherwise ignore
+        if self.grid_jump >= self.num_levels:
+            raise ValueError("grid_jump must be less than num_levels")
+
+        if self.groups is None:
+            sample = self.sample_oat()
+        else:
+            sample = self.sample_groups()
+
+        if self.optimal_trajectories is not None:
             if self.optimal_trajectories >= self.samples:
                 raise ValueError("The number of optimal trajectories should be less than the number of samples.")
             elif self.optimal_trajectories > 10:
                 raise ValueError("Running optimal trajectories greater than values of 10 will take a long time.")
-            elif self.optimal_trajectories <= 1:
+            elif self.optimal_trajectories < 2:
                 raise ValueError("The number of optimal trajectories must be set to 2 or more.")
+            else:
+                sample = self.optimize_trajectories(sample)
 
-        if self.groups is None:
+        self.output_sample = sample
 
-            self.create_sample()
 
-        else:
+    def sample_oat(self):
 
-            self.create_sample_with_groups()
+        D = self.num_vars
+        N = self.samples
 
+        # orientation matrix B: lower triangular (1) + upper triangular (-1)
+        B = np.tril(np.ones([D + 1, D], dtype=int), -1) + \
+            np.triu(-1 * np.ones([D + 1, D], dtype=int))
 
-    def flatten(self, l):
-        for el in l:
-            if isinstance(el, Iterable) and not isinstance(el, str):
-                for sub in self.flatten(el):
-                    yield sub
-            else:
-                yield el
+        # grid step delta, and final sample matrix X
+        delta = self.grid_jump / (self.num_levels - 1)
+        X = np.empty([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)
+        for j in range(N):
 
-    def create_sample(self):
+            # directions matrix DM - diagonal matrix of either +1 or -1
+            DM = np.diag([rd.choice([-1, 1]) for _ in range(D)])
 
-        if self.optimal_trajectories is None:
+            # permutation matrix P
+            perm = np.random.permutation(D)
+            P = np.zeros([D, D])
+            for i in range(D):
+                P[i, perm[i]] = 1
 
-            optimal_sample = morris_oat.sample(self.samples,
-                                               self.parameter_file,
-                                               self.num_levels,
-                                               self.grid_jump)
+            # starting point for this trajectory
+            x_base = np.empty([D + 1, D])
+            for i in range(D):
+                x_base[:, i] = (
+                    rd.choice(np.arange(self.num_levels - self.grid_jump))) / (self.num_levels - 1)
 
-        else:
+            # Indices to be assigned to X, corresponding to this trajectory
+            index_list = np.arange(D + 1) + j * (D + 1)
+            delta_diag = np.diag([delta for _ in range(D)])
 
-            sample = morris_oat.sample(self.samples,
-                                       self.parameter_file,
-                                       self.num_levels,
-                                       self.grid_jump)
-            optimal_sample = \
-                morris_optimal.find_optimum_trajectories(sample,
-                                                         self.samples,
-                                                         self.num_vars,
-                                                         self.optimal_trajectories)
+            X[index_list, :] = 0.5 * \
+                (np.mat(B) * np.mat(P) * np.mat(DM) + 1) * \
+                np.mat(delta_diag) + np.mat(x_base)
 
-        self.output_sample = optimal_sample
+        return X
 
+    def sample_groups(self):
+        '''
+        Returns an N(g+1)-by-k array of N trajectories;
+        where g is the number of groups and k is the number of factors
+        '''
+        N = self.samples
+        G = self.groups
 
-    def create_sample_with_groups(self):
+        if G is None:
+            raise ValueError("Please define the matrix G.")
+        if type(G) is not np.matrixlib.defmatrix.matrix:
+           raise TypeError("Matrix G should be formatted as a numpy matrix")
 
-        self.output_sample = morris_groups.sample(self.samples,
-                                                  self.groups,
-                                                  self.num_levels,
-                                                  self.grid_jump)
-        if self.optimal_trajectories is not None:
-            self.output_sample = \
-                morris_optimal.find_optimum_trajectories(self.output_sample,
-                                                         self.samples,
-                                                         self.num_vars,
-                                                         self.optimal_trajectories)
+        k = G.shape[0]
+        g = G.shape[1]
+        sample = np.empty((N*(g + 1), k))
+        sample = np.array([generate_trajectory(G, self.num_levels, self.grid_jump) for n in range(N)])
+        return sample.reshape((N*(g + 1), k))
+
+
+    def optimize_trajectories(self, input_sample):
+
+        N = self.samples
+        k_choices = self.optimal_trajectories
+
+        if np.any((input_sample < 0) | (input_sample > 1)):
+            raise ValueError("Input sample must be scaled between 0 and 1")
+
+        scores = find_most_distant(input_sample, N, self.num_vars, k_choices)
+
+        index_list = []
+        for j in range(N):
+            index_list.append(np.arange(self.num_vars + 1) + j * (self.num_vars + 1))
+
+        maximum_combo = find_maximum(scores, N, k_choices)
+        output = np.zeros((np.size(maximum_combo) * (self.num_vars + 1), self.num_vars))
+        for counter, x in enumerate(maximum_combo):
+            output[index_list[counter]] = np.array(input_sample[index_list[x]])
+        return output
 
 
     def debug(self):
@@ -138,14 +165,12 @@ def debug(self):
                         default=2, help='Grid jump size (Morris only)')
     parser.add_argument('-k','--k-optimal', type=int, required=False,
                         default=None, help='Number of optimal trajectories (Morris only)')
-    parser.add_argument('--group', type=str, required=False, default=None,
-                       help='File path to grouping file (Morris only)')
     args = parser.parse_args()
 
     np.random.seed(args.seed)
     rd.seed(args.seed)
 
     sample = Morris(args.paramfile, args.samples, args.levels, \
-                    args.grid_jump, args.group, args.k_optimal)
+                    args.grid_jump, args.k_optimal)
 
     sample.save_data(args.output, delimiter=args.delimiter, precision=args.precision)