Source code for SALib.sample.morris.gurobi

'''Finds optimal trajectories using a global optimisation method

Example
-------
Run using

>>> optimal_trajectories.py -n=10 \
    -p=esme_param.txt -o=test_op.txt \
    -s=12892 --num-levels=4 --grid-jump=2 \
    --k-choices=4

'''
from __future__ import division

import re
from datetime import datetime as dt

import numpy as np

try:
    from scipy.misc import comb as nchoosek
except ImportError:
    from scipy.special import comb as nchoosek

from . strategy import Strategy

from SALib.util import requires_gurobipy

try:
    from gurobipy import Model, quicksum, GRB
except ImportError:
    _has_gurobi = False
else:
    _has_gurobi = True


@requires_gurobipy(_has_gurobi)
class GlobalOptimisation(Strategy):
    """Implements the global optimisation algorithm
    """

    def _sample(self, input_sample, num_samples,
                num_params, k_choices, num_groups=None):
        return self.return_max_combo(input_sample, num_samples,
                                     num_params, k_choices, num_groups)

    @staticmethod
    def global_model(N, k_choices, distance_matrix):

        if k_choices >= N:
            raise ValueError("k_choices must be less than N")

        model = Model("distance1")
        trajectories = range(N)

        distance_matrix = np.array(
            distance_matrix / distance_matrix.max(), dtype=np.float64)

        dm = distance_matrix ** 2

        y, x = {}, {}
        for i in trajectories:
            y[i] = model.addVar(vtype="B", obj=0, name="y[%s]" % i)
            for j in range(i + 1, N):
                x[i, j] = model.addVar(
                    vtype="B", obj=1.0, name="x[%s,%s]" % (i, j))
        model.update()

        model.setObjective(quicksum([x[i, j] * dm[j][i]
                                     for i in trajectories
                                     for j in range(i + 1, N)]))

        # Add constraints to the model
        model.addConstr(quicksum([y[i]
                                  for i in trajectories]) <= k_choices, "27")

        for i in trajectories:
            for j in range(i + 1, N):
                model.addConstr(x[i, j] <= y[i], "28-%s-%s" % (i, j))
                model.addConstr(x[i, j] <= y[j], "29-%s-%s" % (i, j))
                model.addConstr(y[i] + y[j] <= 1 + x[i, j],
                                "30-%s-%s" % (i, j))

        model.addConstr(quicksum([x[i, j] for i in trajectories
                                  for j in range(i + 1, N)])
                        <= nchoosek(k_choices, 2), "Cut_1")
        model.update()
        return model

    def return_max_combo(self, input_data, N, num_params,
                         k_choices, num_groups=None):
        """Find the optimal combination of most different trajectories

        Arguments
        ---------
        input_data : numpy.ndarray
        N : int
            The number of samples
        num_params : int
            The number of factors
        k_choices : int
            The number of optimal trajectories to select
        num_groups: int, default=None
            The number of groups

        Returns
        -------
        list
        """

        distance_matrix = self.compute_distance_matrix(
            input_data, N, num_params, num_groups)

        model = self.global_model(N, k_choices, distance_matrix)
        model.params.Threads = 0

        model.ModelSense = GRB.MAXIMIZE
        model.optimize()
        if model.Status == GRB.OPTIMAL:
            print(model.objval)

        variables = list(model.getVars())
        x_vars = []
        for v in variables:
            if (v.X > 0) & (v.VarName[0] == 'x'):
                x_vars.append(v.VarName)
        b = [re.findall(r"\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)


[docs]def timestamp(num_params, p_levels, k_choices, N): """ Returns a uniform timestamp with parameter values for file identification """ string = "_v%s_l%s_gs%s_k%s_N%s_%s.txt" % (num_params, p_levels, k_choices, N, dt.strftime(dt.now(), "%d%m%y%H%M%S")) return string