text: """
  Generate a sample using the Method of Morris
  
  Three variants of Morris' sampling for elementary effects is supported:
  
  - Vanilla Morris
  - Optimised trajectories when ``optimal_trajectories=True`` (using
      Campolongo's enhancements from 2007 and optionally Ruano's enhancement
      from 2012; ``local_optimization=True``)
  - Groups with optimised trajectories when ``optimal_trajectories=True`` and
      the problem definition specifies groups (note that ``local_optimization``
      must be ``False``)
  
  At present, optimised trajectories is implemented using either a brute-force
  approach, which can be very slow, especially if you require more than four
  trajectories, or a local method based which is much faster. Both methods now
  implement working with groups of factors.
  
  Note that the number of factors makes little difference,
  but the ratio between number of optimal trajectories and the sample size
  results in an exponentially increasing number of scores that must be
  computed to find the optimal combination of trajectories.  We suggest going
  no higher than 4 from a pool of 100 samples with the brute force approach.
  With local_optimization = True (which is default),
  it is possible to go higher than the previously suggested 4 from 100.
  
  """
  from __future__ import division
  
  import numpy as np
  
  import numpy.random as rd
  import warnings
  
  from . gurobi import GlobalOptimisation
  from . local import LocalOptimisation
  from . brute import BruteForce
  
  from . strategy import SampleMorris
  
  from SALib.sample import common_args
  from SALib.util import scale_samples, read_param_file, compute_groups_matrix
  
  try:
      import gurobipy
  except ImportError:
      _has_gurobi = False
  else:
      _has_gurobi = True
  
  __all__ = ['sample']
  
  
  def sample(problem, N, num_levels=4, optimal_trajectories=None,
             local_optimization=True, seed=None):
      """Generate model inputs using the Method of Morris
  
      Returns a NumPy matrix containing the model inputs required for Method of
      Morris.  The resulting matrix has :math:`(G+1)*T` rows and :math:`D`
      columns, where :math:`D` is the number of parameters, :math:`G` is the
      number of groups (if no groups are selected, the number of parameters).
      :math:`T` is the number of trajectories :math:`N`,
      or `optimal_trajectories` if selected.
      These model inputs  are intended to be used with
      :func:`SALib.analyze.morris.analyze`.
  
      Parameters
      ----------
      problem : dict
          The problem definition
      N : int
          The number of trajectories to generate
      num_levels : int, default=4
          The number of grid levels (should be even)
      optimal_trajectories : int
          The number of optimal trajectories to sample (between 2 and N)
      local_optimization : bool, default=True
          Flag whether to use local optimization according to Ruano et al. (2012)
          Speeds up the process tremendously for bigger N and num_levels.
          If set to ``False`` brute force method is used, unless ``gurobipy`` is
          available
  
      Returns
      -------
      sample : numpy.ndarray
          Returns a numpy.ndarray containing the model inputs required for Method
          of Morris. The resulting matrix has :math:`(G/D+1)*N/T` rows and
          :math:`D` columns, where :math:`D` is the number of parameters.
      """
      if seed:
          np.random.seed(seed)
  
      if not num_levels % 2 == 0:
          warnings.warn("num_levels should be an even number, sample may be biased")
      if problem.get('groups'):
          sample = _sample_groups(problem, N, num_levels)
      else:
          sample = _sample_oat(problem, N, num_levels)
  
      if optimal_trajectories:
  
          sample = _compute_optimised_trajectories(problem,
                                                   sample,
                                                   N,
                                                   optimal_trajectories,
                                                   local_optimization)
  
      scale_samples(sample, problem['bounds'])
      return sample
  
  
  def _sample_oat(problem, N, num_levels=4):
      """Generate trajectories without groups
  
      Arguments
      ---------
      problem : dict
          The problem definition
      N : int
          The number of samples to generate
      num_levels : int, default=4
          The number of grid levels
      """
      group_membership = np.asmatrix(np.identity(problem['num_vars'],
                                                 dtype=int))
  
      num_params = group_membership.shape[0]
      sample = np.array([generate_trajectory(group_membership,
                                             num_levels)
                         for n in range(N)])
      return sample.reshape((N * (num_params + 1), num_params))
  
  
  def _sample_groups(problem, N, num_levels=4):
      """Generate trajectories for groups
  
      Returns an :math:`N(g+1)`-by-:math:`k` array of `N` trajectories,
      where :math:`g` is the number of groups and :math:`k` is the number
      of factors
  
      Arguments
      ---------
      problem : dict
          The problem definition
      N : int
          The number of trajectories to generate
      num_levels : int, default=4
          The number of grid levels
  
      Returns
      -------
      numpy.ndarray
      """
      if len(problem['groups']) != problem['num_vars']:
          raise ValueError("Groups do not match to number of variables")
  
      group_membership, _ = compute_groups_matrix(problem['groups'])
  
      if group_membership is None:
          raise ValueError("Please define the 'group_membership' matrix")
      if not isinstance(group_membership, np.ndarray):
          raise TypeError("Argument 'group_membership' should be formatted \
                           as a numpy ndarray")
  
      num_params = group_membership.shape[0]
      num_groups = group_membership.shape[1]
      sample = np.zeros((N * (num_groups + 1), num_params))
      sample = np.array([generate_trajectory(group_membership,
                                             num_levels)
                         for n in range(N)])
      return sample.reshape((N * (num_groups + 1), num_params))
  
  
  def generate_trajectory(group_membership, num_levels=4):
      """Return a single trajectory
  
      Return a single trajectory of size :math:`(g+1)`-by-:math:`k`
      where :math:`g` is the number of groups,
      and :math:`k` is the number of factors,
      both implied by the dimensions of `group_membership`
  
      Arguments
      ---------
      group_membership : np.ndarray
          a k-by-g matrix which notes factor membership of groups
      num_levels : int, default=4
          The number of levels in the grid
  
      Returns
      -------
      np.ndarray
      """
  
      delta = compute_delta(num_levels)
  
      # Infer number of groups `g` and number of params `k` from
      # `group_membership` matrix
      num_params = group_membership.shape[0]
      num_groups = group_membership.shape[1]
  
      # Matrix B - size (g + 1) * g -  lower triangular matrix
      B = np.tril(np.ones([num_groups + 1, num_groups],
                          dtype=int), -1)
  
      P_star = generate_p_star(num_groups)
  
      # Matrix J - a (g+1)-by-num_params matrix of ones
      J = np.ones((num_groups + 1, num_params))
  
      # Matrix D* - num_params-by-num_params matrix which decribes whether
      # factors move up or down
      D_star = np.diag(rd.choice([-1, 1], num_params))
  
      x_star = generate_x_star(num_params, num_levels)
  
      # Matrix B* - size (num_groups + 1) * num_params
      B_star = compute_b_star(J, x_star, delta, B,
                              group_membership, P_star, D_star)
  
      return B_star
  
  
  def compute_b_star(J, x_star, delta, B, G, P_star, D_star):
      """
      """
      element_a = J[0, :] * x_star
      element_b = np.matmul(G, P_star).T
      element_c = np.matmul(2.0 * B, element_b)
      element_d = np.matmul((element_c - J), D_star)
  
      b_star = element_a + (delta / 2.0) * (element_d + J)
      return b_star
  
  
  def generate_p_star(num_groups):
      """Describe the order in which groups move
  
      Arguments
      ---------
      num_groups : int
  
      Returns
      -------
      np.ndarray
          Matrix P* - size (g-by-g)
      """
      p_star = np.eye(num_groups, num_groups)
      rd.shuffle(p_star)
      return p_star
  
  
  def generate_x_star(num_params, num_levels):
      """Generate an 1-by-num_params array to represent initial position for EE
  
      This should be a randomly generated array in the p level grid
      :math:`\omega`
  
      Arguments
      ---------
      num_params : int
          The number of parameters (factors)
      num_levels : int
          The number of levels
  
      Returns
      -------
      numpy.ndarray
          The initial starting positions of the trajectory
      """
      x_star = np.zeros((1, num_params))
      delta = compute_delta(num_levels)
      bound = 1 - delta
      grid = np.linspace(0, bound, int(num_levels / 2))
  
      x_star[0, :] = rd.choice(grid, num_params)
  
      return x_star
  
  
  def compute_delta(num_levels):
      """Computes the delta value from number of levels
  
      Arguments
      ---------
      num_levels : int
          The number of levels
  
      Returns
      -------
      float
      """
      return num_levels / (2.0 * (num_levels - 1))
  
  
  def _compute_optimised_trajectories(problem, input_sample, N, k_choices,
                                      local_optimization=False):
      '''
      Calls the procedure to compute the optimum k_choices of trajectories
      from the input_sample.
      If there are groups, then this procedure allocates the groups to the
      correct call here.
  
      Arguments
      ---------
      problem : dict
          The problem definition
      input_sample :
      N : int
          The number of samples to generate
      k_choices : int
          The number of optimal trajectories
      local_optimization : bool, default=False
          If true, uses local optimisation heuristic
      '''
      if _has_gurobi is False \
              and local_optimization is False \
              and k_choices > 10:
          msg = "Running optimal trajectories greater than values of 10 \
                  will take a long time."
          raise ValueError(msg)
  
      num_params = problem['num_vars']
  
      if np.any((input_sample < 0) | (input_sample > 1)):
          raise ValueError("Input sample must be scaled between 0 and 1")
  
      if _has_gurobi and local_optimization is False:
          # Use global optimization method
          strategy = GlobalOptimisation()
      elif local_optimization:
          # Use local method
          strategy = LocalOptimisation()
      else:
          # Use brute force approach
          strategy = BruteForce()
  
      if problem.get('groups'):
          num_groups = len(set(problem['groups']))
      else:
          num_groups = None
  
      context = SampleMorris(strategy)
      output = context.sample(input_sample, N, num_params,
                              k_choices, num_groups)
  
      return output
  
  
  def cli_parse(parser):
      parser.add_argument('-l', '--levels', type=int, required=False,
                          default=4, help='Number of grid levels \
                          (Morris only)')
      parser.add_argument('-k', '--k-optimal', type=int, required=False,
                          default=None,
                          help='Number of optimal trajectories \
                          (Morris only)')
      parser.add_argument('-lo', '--local', type=bool, required=True,
                          default=False,
                          help='Use the local optimisation method \
                          (Morris with optimization only)')
      return parser
  
  
  def cli_action(args):
      rd.seed(args.seed)
  
      problem = read_param_file(args.paramfile)
      param_values = sample(problem, args.samples, args.levels,
                            args.k_optimal, args.local)
  
      np.savetxt(args.output, param_values, delimiter=args.delimiter,
                 fmt='%.' + str(args.precision) + 'e')
  
  
  if __name__ == "__main__":
      common_args.run_cli(cli_parse, cli_action)
  