def K_dominnace_check_2(self, u_d, v_d, _inequalities):
        """

        :param u_d: a d-dimensional vector(list) like [ 8.53149891  3.36436796]
        :param v_d: tha same list like u_d
        :param _inequalities: list of constraints on d-dimensional Lambda Polytope like
         [[0, 1, 0], [1, -1, 0], [0, 0, 1], [1, 0, -1], [0.0, 1.4770889, -3.1250839]]
        :return: True if u is Kdominance to v regarding given _inequalities otherwise False
        """
        _d = len(u_d)

        prob = LpProblem("Kdominance", LpMinimize)
        lambda_variables = LpVariable.dicts("l", range(_d), 0)

        for inequ in _inequalities:
            prob += lpSum([inequ[j + 1] * lambda_variables[j] for j in range(0, _d)]) + inequ[0] >= 0

        prob += lpSum([lambda_variables[i] * (u_d[i]-v_d[i]) for i in range(_d)])

        #prob.writeLP("show-Ldominance.lp")

        status = prob.solve()
        LpStatus[status]

        result = value(prob.objective)
        if result < 0:
            return False

        return True

def solve(g):
    el = g.get_edge_list()
    nl = g.get_node_list()
    p = LpProblem('min_cost', LpMinimize)
    capacity = {}
    cost = {}
    demand = {}
    x = {}
    for e in el:
        capacity[e] = g.get_edge_attr(e[0], e[1], 'capacity')
        cost[e] = g.get_edge_attr(e[0], e[1], 'cost')
    for i in nl:
        demand[i] = g.get_node_attr(i, 'demand')
    for e in el:
        x[e] = LpVariable("x"+str(e), 0, capacity[e])
    # add obj
    objective = lpSum (cost[e]*x[e] for e in el)
    p += objective
    # add constraints
    for i in nl:
        out_neig = g.get_out_neighbors(i)
        in_neig = g.get_in_neighbors(i)
        p += lpSum(x[(i,j)] for j in out_neig) -\
             lpSum(x[(j,i)] for j in in_neig)==demand[i]
    p.solve()
    return x, value(objective)

def pi_solve(pi, Q, s):
    # The 'prob' variable will contain the problem data.
    prob = LpProblem('FindPi', LpMinimize)

    a0 = LpVariable('a0', 0.0)  # the minimum is 0.0
    a1 = LpVariable('a1', 0.0)  # the minimum is 0.0
    a2 = LpVariable('a2', 0.0)  # the minimum is 0.0
    a3 = LpVariable('a3', 0.0)  # the minimum is 0.0
    a4 = LpVariable('a4', 0.0)  # the minimum is 0.0
    v = LpVariable('v', 0.0)

    # The objective function is added to 'prob' first
    prob += v, "to minimize"

    # constraints
    prob +=  a0 * Q[s,0,0] + a1 * Q[s,1,0] + a2 * Q[s,2,0] + a3 * Q[s,3,0] + a4 * Q[s,4,0] <= v, 'constraint 1'
    prob +=  a0 * Q[s,0,1] + a1 * Q[s,1,1] + a2 * Q[s,2,1] + a3 * Q[s,3,1] + a4 * Q[s,4,1] <= v, 'constraint 2'
    prob +=  a0 * Q[s,0,2] + a1 * Q[s,1,2] + a2 * Q[s,2,2] + a3 * Q[s,3,2] + a4 * Q[s,4,2] <= v, 'constraint 3'
    prob +=  a0 * Q[s,0,3] + a1 * Q[s,1,3] + a2 * Q[s,2,3] + a3 * Q[s,3,3] + a4 * Q[s,4,3] <= v, 'constraint 4'
    prob +=  a0 * Q[s,0,4] + a1 * Q[s,1,4] + a2 * Q[s,2,4] + a3 * Q[s,3,4] + a4 * Q[s,4,4] <= v, 'constraint 5'

    prob +=  1.0*a0 + 1.0*a1 + 1.0*a2 + 1.0*a3 + 1.0*a4 == 1, 'constraint 6'

    prob.solve()

    pi_prime = [a.varValue for a in prob.variables()[:5]]

    return np.array(pi_prime)

    def K_dominance_check(self, _V_best_d, Q_d):
        """
        :param _V_best_d: a list of d-dimension
        :param Q_d: a list of d-dimension
        :return: True if _V_best_d is prefered to Q_d regarding self.Lambda_inequalities and using Kdominance
         other wise it returns False
        """
        _d = len(_V_best_d)

        prob = LpProblem("Ldominance", LpMinimize)
        lambda_variables = LpVariable.dicts("l", range(_d), 0)

        for inequ in self.Lambda_ineqalities:
            prob += lpSum([inequ[j + 1] * lambda_variables[j] for j in range(0, _d)]) + inequ[0] >= 0

        prob += lpSum([lambda_variables[i] * (_V_best_d[i]-Q_d[i]) for i in range(_d)])

        #prob.writeLP("show-Ldominance.lp")

        status = prob.solve()
        LpStatus[status]
        result = value(prob.objective)

        if result < 0:
            return False

        return True

def make_into_lp_problem(good_for, N):
    """
    Helper function for solve_with_lp_and_reduce()

    N --- number of isoform sequences
    good_for --- dict of <isoform_index> --> list of matched paths index
    """
    prob = LpProblem("The Whiskas Problem",LpMinimize)

    # each good_for is (isoform_index, [list of matched paths index])
    # ex: (0, [1,2,4])
    # ex: (3, [2,5])
    variables = [LpVariable(str(i),0,1,LpInteger) for i in xrange(N)]

    # objective is to minimize sum_{Xi}
    prob += sum(v for v in variables)

    # constraints are for each isoform, expressed as c_i * x_i >= 1
    # where c_i = 1 if x_i is matched for the isoform
    # ex: (0, [1,2,4]) becomes t_0 = x_1 + x_2 + x_4 >= 1
    for t_i, p_i_s in good_for:
        #c_i_s = [1 if i in p_i_s else 0 for i in xrange(N)]
        prob += sum(variables[i]*(1 if i in p_i_s else 0) for i in xrange(N)) >= 1
    prob.writeLP('cogent.lp')
    return prob

	def get_linear_program_solution(self, c, b, A, x):
		prob = LpProblem("myProblem", LpMinimize)
		prob += lpSum([xp*cp for xp, cp in zip(x, c)]), "Total Cost of Ingredients per can" 

		for row, cell in zip(A, b):
			prob += lpSum(ap*xp for ap, xp in zip(row,  x)) <= cell

		solved = prob.solve()
		return prob

def fuelSolution(airports):
    for i in range(len(airports)):
        airports[i].previous = airports[i-1]
    problem = LpProblem('Flight', LpMinimize)
    problem+=sum(airport.fuelcost*airport.refuel for airport in airports), 'Cost'
    for ap in airports:
        problem+= ap.min_fuel + ap.refuel <= ap.startingFuel, 'Min Fuel %s' % ap.startingFuel
        problem+= ap.startingFuel == ap.previous.startingFuel + ap.refuel - ap.distance*(1 + (ap.startingFuel + ap.previous.startingFuel - ap.refuel)/(2.*rate)), 'Takeoff Fuel Level %s' % ap.startingFuel
    problem.solve()
    return problem

def main():
    x = LpVariable("x", 0, 3)
    y = LpVariable("y", 0, 1)
    prob = LpProblem("myProblem", LpMinimize)
    prob += x + y <= 2
    prob += -4*x + y
    status = prob.solve(COIN(msg = 0))
    print LpStatus[status]
    print value(x)
    print value(y)

def knapsack01(obj, weights, capacity):
    """ 0/1 knapsack solver, maximizes profit. weights and capacity integer """
        
    debug_subproblem = False
    
    assert len(obj) == len(weights)
    n = len(obj)
    if n == 0:
        return 0, []

    if debug_subproblem:
        relaxation = LpProblem('relaxation', LpMaximize)
        relax_vars = [str(i) for i in range(n)]
        var_dict   = LpVariable.dicts("", relax_vars, 0, 1, LpBinary)
        relaxation += (lpSum(var_dict[str(i)] * weights[i] for i in range(n)) 
                       <= capacity)
        relaxation += lpSum(var_dict[str(i)] * obj[i] for i in range(n))
        relaxation.solve()
        relax_obj = value(relaxation.objective)

        solution =  [i for i in range(n) if var_dict[str(i)].varValue > tol ]

        print relax_obj, solution


    c = [[0]*(capacity+1) for i in range(n)]
    added = [[False]*(capacity+1) for i in range(n)]
    # c [items, remaining capacity]
    # important: this code assumes strictly positive objective values
    for i in range(n):
        for j in range(capacity+1):
            if (weights[i] > j):
                c[i][j] = c[i-1][j]
            else:
                c_add = obj[i] + c[i-1][j-weights[i]]
                if c_add > c[i-1][j]:
                    c[i][j] = c_add
                    added[i][j] = True
                else:
                    c[i][j] = c[i-1][j]

    # backtrack to find solution
    i = n-1
    j = capacity

    solution = []
    while i >= 0 and j >= 0:
        if added[i][j]:
            solution.append(i)
            j -= weights[i]
        i -= 1
        
    return c[n-1][capacity], solution

def pe185():
    """
    Modelling as an integer programming problem.
    Then using PuLP to solve it. It's really fast, just 0.24 seconds. 
    For details, see https://pythonhosted.org/PuLP/index.html
    """
    
    from pulp import LpProblem, LpVariable, LpMinimize, LpInteger, lpSum, value

    constraints = [
        ('2321386104303845', 0),
        ('3847439647293047', 1),
        ('3174248439465858', 1),
        ('8157356344118483', 1),
        ('6375711915077050', 1),
        ('6913859173121360', 1),
        ('4895722652190306', 1),
        ('5616185650518293', 2),
        ('4513559094146117', 2),
        ('2615250744386899', 2),
        ('6442889055042768', 2),
        ('2326509471271448', 2),
        ('5251583379644322', 2),
        ('2659862637316867', 2),
        ('5855462940810587', 3),
        ('9742855507068353', 3),
        ('4296849643607543', 3),
        ('7890971548908067', 3),
        ('8690095851526254', 3),
        ('1748270476758276', 3),
        ('3041631117224635', 3),
        ('1841236454324589', 3)
    ]

    VALs = map(str, range(10))
    LOCs = map(str, range(16))
    choices = LpVariable.dicts("Choice", (LOCs, VALs), 0, 1, LpInteger)

    prob = LpProblem("pe185", LpMinimize)
    prob += 0, "Arbitrary Objective Function"

    for s in LOCs:
        prob += lpSum([choices[s][v] for v in VALs]) == 1, ""

    for c, n in constraints:
        prob += lpSum([choices[str(i)][v] for i,v in enumerate(c)]) == n, ""

    prob.writeLP("pe185.lp")
    prob.solve()
    res = int(''.join(v for s in LOCs for v in VALs if value(choices[s][v])))

    # answer: 4640261571849533
    return res

 def get_best_assignments(self, row):
     df = self.assignments[(self.assignments.FromIcao == row['FromIcao']) &
                           (self.assignments.ToIcao == row['ToIcao']) & (self.assignments.Amount <= row['Seats'])]
     if not len(df):
         return None
     prob = LpProblem("Knapsack problem", LpMaximize)
     w_list = df.Amount.tolist()
     p_list = df.Pay.tolist()
     x_list = [LpVariable('x{}'.format(i), 0, 1, 'Integer') for i in range(1, 1 + len(w_list))]
     prob += sum([x * p for x, p in zip(x_list, p_list)]), 'obj'
     prob += sum([x * w for x, w in zip(x_list, w_list)]) <= row['Seats'], 'c1'
     prob.solve()
     return df.iloc[[i for i in range(len(x_list)) if x_list[i].varValue]]

def solve_under_coverage(graph, min_coverage=80):

    prob = LpProblem("granularity selection", LpMinimize)
    codelet_vars = LpVariable.dicts("codelet",
            graph,
            lowBound=0,
            upBound=1,
            cat=LpInteger)

    # Objective function: minimize the total replay cost of selected codelets

    # Compute replay time
    for n,d in graph.nodes(data=True):
      d['_total_replay_cycles'] = 0
      for inv in d['_invocations']:
        d['_total_replay_cycles'] = d['_total_replay_cycles'] + float(inv["Invivo (cycles)"])

    prob += lpSum([codelet_vars[n]*d['_total_replay_cycles'] for n,d in graph.nodes(data=True)])

    # and with good coverage
    prob += (lpSum([codelet_vars[n]*d['_coverage'] for n,d in graph.nodes(data=True)]) >= min_coverage)

    # selected codelets should match
    for n,d in graph.nodes(data=True):
        if not d['_matching']:
            prob += codelet_vars[n] == 0

    # Finally we should never include both the children and the parents
    for dad in graph.nodes():
        for son in graph.nodes():
            if not dad in nx.ancestors(graph, son):
                continue
            # We cannot select dad and son at the same time
            prob += codelet_vars[dad] + codelet_vars[son] <= 1

    #prob.solve(GLPK())
    prob.solve()
    if (LpStatus[prob.status] != 'Optimal'):
        raise Unsolvable()

    for v in prob.variables():
        assert v.varValue == 1.0 or v.varValue == 0.0
        if v.varValue == 1.0:

            for n,d in graph.nodes(data=True):
                if ("codelet_"+str(n)) == v.name:
                    d["_selected"] = True
                    yield n

 def setup(self, statement):
     self.problem = LpProblem(statement, LpMinimize)
     self.problem += sum(self.variables.values())
     
     for a, b in combinations(self.candidates, 2):
         ab = a + b
         ba = b + a
         
         if ab in statement and ba in statement:
             # The statement itself contains the ordering information.
             pass
         elif ab in statement:
             self.constrainGreater(ab, ba)
         elif ba in statement:
             self.constrainGreater(ba, ab)
         #else:
         #    self.constrainEqual(ab, ba)
     
     prev = None
     for rank in statement.split(">"):
         pairs = rank.split("=")
         if prev:
             self.constrainGreater(prev, pairs[0])
         for n in range(len(pairs) - 1):
             self.constrainEqual(pairs[n], pairs[n+1])
         prev = pairs[-1]

def get_optimal_routes(sources, destinations):
    sources = collections.OrderedDict([(x['id'], x) for x in sources])
    destinations = collections.OrderedDict([(x['id'], x) for x in destinations])

    sources_points = [{'lat': x['lat'], 'lng': x['lng']} for x in sources.values()]
    destinations_points = [{'lat': x['lat'], 'lng': x['lng']} for x in destinations.values()]

    source_ids = [str(x['id']) for x in sources.values()]
    dest_ids = [str(x['id']) for x in destinations.values()]

    demand = {str(x['id']): convert_int(x['num_students']) for x in sources.values()}
    supply = {str(x['id']): convert_int(x['num_students']) for x in destinations.values()}

    log.info("Calling gmaps api...")
    distances = gmaps.distance_matrix(origins=sources_points, destinations=destinations_points, mode='walking')

    costs = {}
    for i, origin in enumerate(distances['rows']):
        origin_costs = {}
        for j, entry in enumerate(origin['elements']):
            origin_costs[dest_ids[j]] = entry['duration']['value']
        costs[source_ids[i]] = origin_costs

    prob = LpProblem("Evaucation Routing for Schools",LpMinimize)
    routes = [(s,d) for s in source_ids for d in dest_ids]
    route_lookup = {'Route_{}_{}'.format(x.replace(' ','_'),y.replace(' ','_')):(x,y) for (x,y) in routes}
    route_vars = LpVariable.dicts("Route",(source_ids,dest_ids),0,None,LpInteger)
    prob += lpSum([route_vars[w][b]*(costs[w][b]**2) for (w,b) in routes])
    for dest in dest_ids:
        prob += lpSum([route_vars[source][dest] for source in source_ids]) <= supply[dest], "Students going to {} is <= {}".format(dest, supply[dest])
    for source in source_ids:
        prob += lpSum([route_vars[source][dest] for dest in dest_ids]) == demand[source], "Students leaving {} is {}".format(source, demand[source])

    log.info("Optimizing routes...")
    prob.solve()

    if prob.status != 1:
        raise Exception("Algorithm could not converge to a solution")

    result = []
    for v in prob.variables():
        src, dst = route_lookup[v.name]
        value = v.value()
        result.append({'src': sources[src], 'dst': destinations[dst], 'value': int(value)})
    return result

    def _recurse(work_pairs, mae, grace=None):
        req = grace or 1.0

        demands = {
            'morning': mae[0],
            'afternoon': mae[1],
            'exam': mae[2],
        }
        total_avg = float(sum([demands[i] for i in SHIFTS])) / work_pairs
        total_low, total_high = floor(total_avg), ceil(total_avg)
        work_pair_count = work_pairs
        avgs = [float(demands[i]) / work_pair_count for i in SHIFTS]
        lows = [floor(a) for a in avgs]
        highs = [ceil(a) for a in avgs]

        target = req * total_avg * float(sum([COSTS[i] for i in SHIFTS])) / len(SHIFTS)

        prob = LpProblem("Work Distribution", LpMinimize)
        var_prefix = "shift"
        shift_vars = LpVariable.dicts(var_prefix, SHIFTS, 0, cat=LpInteger)
        prob += lpSum([COSTS[i] * shift_vars[i] for i in SHIFTS]), "cost of combination"
        prob += lpSum([COSTS[i] * shift_vars[i] for i in SHIFTS]) >= target, "not too good"
        prob += lpSum([shift_vars[i] for i in SHIFTS]) >= total_low, "low TOTAL"
        prob += lpSum([shift_vars[i] for i in SHIFTS]) <= total_high, "high TOTAL"

        for shift, low, high in zip(SHIFTS, lows, highs):
            prob += lpSum([shift_vars[shift]]) >= low, "low %s" % shift
            prob += lpSum([shift_vars[shift]]) <= high, "high %s" % shift

        prob.solve(GLPK_CMD(msg=0))

        if not LpStatus[prob.status] == 'Optimal':
            next_grace = req - 0.1
            assert 0.0 < next_grace
            return _recurse(work_pairs, mae, next_grace)

        new_mae = [0, 0, 0]
        solution = [0, 0, 0]
        for v in prob.variables():
            for pos, name in enumerate(SHIFTS):
                if v.name == "%s_%s" % (var_prefix, name):
                    solution[pos] = v.varValue
                    new_mae[pos] = mae[pos] - solution[pos]

        return (PairAlloc(solution), work_pairs - 1) + tuple(new_mae)

 def create_prob(self, sense=LpMaximize, use_glpk=False):
     # create the LP
     self.prob = LpProblem('OptKnock', sense=sense)
     if use_glpk:
         self.prob.solver = solvers.GLPK()
     else:
         self.prob.solver = solvers.CPLEX(msg=self.verbose)
     if not self.prob.solver.available():
         raise Exception("CPLEX not available")    

def runLPP(lines):
    # reads the data from the CSV file into the LPP
    variables, ingredients, servingsPerBlock, costsPerServing = [], [], [], []
    
    for row in lines:  # read in the variables from the lines
        variables.append(LpVariable("x_" + row[1], 0, None, LpInteger))
        ingredients.append(row[0])
        servingsPerBlock.append(float(row[2]))
        costsPerServing.append(float(row[3]))
    variables.append(LpVariable("s", 0))

    # read in the upper and lower bounds for aLCM
    minServe = lines[0][4]
    maxServe = lines[0][5]

    # makes the new LP Problem
    problem = LpProblem("Approximate LCM", LpMinimize)

    # serving constraints: specifies the interval in which the approximate
    # LCM may lie
    min_constraint = variables[-1] >= float(minServe)
    max_constraint = variables[-1] <= float(maxServe)
    problem += min_constraint
    problem += max_constraint

    # block constraints: there must be enough of each ingredient for the
    # ultimate optimal number
    for v in range(len(variables) - 1):
        min_constraint = variables[v] * servingsPerBlock[v] >= variables[-1]
        problem += min_constraint

    pps = 0.0  # price per serving
    for i in range(len(variables) - 1):
        pps += costsPerServing[i]

    #  add objective function to be minimized (waste)
    obj_fn = 0 * variables[0]
    for i in range(len(variables) - 1):
        obj_fn += variables[i] * costsPerServing[i] * servingsPerBlock[i]
    obj_fn -= variables[-1] * pps
    problem += obj_fn  # add the objective funciton to the problem
    problem.solve()  # runs the linear programming problem

    return [variables, costsPerServing, obj_fn]

def find_optimal_strategy(states, controls, costs, kernels, solver=None):
    """
    :param states: Number of system states (X)
    :param controls: Number of system controls (U)
    :param сosts: Cost matrix |X| x |U|
    :param kernels: Transition kernels. Dimensionality |X| x |X| x |U|
    """
    tolerance = 10e-15

    X = range(states)
    U = range(controls)
    R = costs
    Q = kernels

    # Check costs
    # Check num of rows
    assert(len(R) == states)
    for row in R:
        # Check num of cols
        assert(len(row) == controls)

    # Check kernels
    # Check num of rows
    assert(len(Q) == states)
    for row in Q:
        # Check num of cols
        assert(len(row) == states)
        for items in row:
            # Check num of items
            assert(len(items) == controls)
        # Check if distribution is normed
        for dist in zip(*row):
            assert(sum(dist)-1 < tolerance)

    # LP object
    optm = LpProblem("Optimal strategy", sense=LpMinimize)

    # Variables (continuous in range [0, 1])
    Z = [[LpVariable("z({},{})".format(x, u), 0, 1) \
                for u in U] for x in X]

    # Objective
    optm.objective = sum(np.dot(Z[x], R[x]) for x in X)

    # Constraints
    for x in X:
        cn = (sum(Z[x]) == sum(Q[y][x][u]*Z[y][u] for u in U for y in X))
        optm.add(cn)
    cn = sum(Z[x][u] for u in U for x in X) == 1
    optm.add(cn)

    optm.solve(solver)

    return [(x, u) for u in U for x in X if value(Z[x][u]) != 0]

	def lp(self): 
		a_name = []
		p_name = []
		n_name = []
		a = []
		p = []
		n = []
		for i in range(self.att):
			a_name.append('a_' + str(i))
		for h in range(self.h):
			p_name.append('p_' + str(h))
		for m in range(self.m):
			n_name.append('n_' + str(m))

		a = [LpVariable(a_name[i]) for i in range(self.att)]
		b = LpVariable('b')

		p = [LpVariable(p_name[i],0) for i in range(self.h)]
		n = [LpVariable(n_name[i],0) for i in range(self.m)]

		#set objective function and constrains
		prob = LpProblem("myProblem",LpMinimize)
		sumP = 0
		sumN = 0
		for i in range(self.h):
			sumX = 0
			for j in range(self.att):
				sumX -= a[j] * self.H[i,j]
			prob += sumX + b + 1 <= p[i]
			sumP += p[i]
		for i in range(self.m):
			sumX = 0
			for j in range(self.att):
				sumX += a[j] * self.M[i,j]
			prob += sumX - b + 1 <= n[i]
			sumN += n[i]
	
		prob += (1/h) * sumP + (1/m) * sumN
		status = prob.solve()

		for i in range(self.att):
			self.a.append(pulp.value(a[i]))
		self.b = pulp.value(b)

def pulp_solve(A, b, c):
    problem = LpProblem(sense = pulp.constants.LpMaximize)

    x = [LpVariable('x' + str(i + 1),
                    0,
                    None,
                    LpInteger)
         for i in xrange(len(c))]

    problem += lpSum(ci * xi for ci, xi in zip(c, x))

    for ai, bi in zip(A, b):
        problem += lpSum(aij * xj for aij, xj in zip(ai, x)) == bi
    
    status = problem.solve()
    
    return (pulp.constants.LpStatus[status],
            [variable.varValue for variable in problem.sortedVariables()],
            problem.objective.value())