def gcp_fixed_k(V,E,K):
    """gcp_fixed_k -- model for minimizing number of bad edges in coloring a graph
    Parameters:
        - V: set/list of nodes in the graph
        - E: set/list of edges in the graph
        - K: number of colors to be used
    Returns a model, ready to be solved.
    """
    model = Model("gcp - fixed k")

    x,z = {},{}
    for i in V:
        for k in range(K):
            x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k))
    for (i,j) in E:
        z[i,j] = model.addVar(vtype="B", name="z(%s,%s)"%(i,j))

    for i in V:
        model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)" % i)

    for (i,j) in E:
        for k in range(K):
            model.addCons(x[i,k] + x[j,k] <= 1 + z[i,j], "BadEdge(%s,%s,%s)"%(i,j,k))

    model.setObjective(quicksum(z[i,j] for (i,j) in E), "minimize")

    model.data = x,z
    return model

def sils(T,f,c,d,h):
    """sils -- LP lotsizing for the single item lot sizing problem
    Parameters:
        - T: number of periods
        - P: set of products
        - f[t]: set-up costs (on period t)
        - c[t]: variable costs
        - d[t]: demand values
        - h[t]: holding costs
    Returns a model, ready to be solved.
    """
    model = Model("single item lotsizing")
    Ts = range(1,T+1)
    M = sum(d[t] for t in Ts)
    y,x,I = {},{},{}
    for t in Ts:
        y[t] = model.addVar(vtype="I", ub=1, name="y(%s)"%t)
        x[t] = model.addVar(vtype="C", ub=M, name="x(%s)"%t)
        I[t] = model.addVar(vtype="C", name="I(%s)"%t)
    I[0] = 0

    for t in Ts:
        model.addCons(x[t] <= M*y[t], "ConstrUB(%s)"%t)
        model.addCons(I[t-1] + x[t] == I[t] + d[t], "FlowCons(%s)"%t)

    model.setObjective(\
        quicksum(f[t]*y[t] + c[t]*x[t] + h[t]*I[t] for t in Ts),\
        "minimize")

    model.data = y,x,I
    return model

def _init():
    model = Model()
    model.hideOutput()
    x = model.addVar("x","B")
    y = model.addVar("y","B")
    z = model.addVar("z","B")
    return model, x, y, z

def create_tsp(vertices, distance):
    model = Model("TSP")

    x = {}  # binary variable to select edges
    for (i, j) in pairs(vertices):
        x[i, j] = model.addVar(vtype = "B", name = "x(%s,%s)" % (i, j))

    for i in vertices:
        model.addCons(quicksum(x[j, i] for j in vertices if j < i) +
                      quicksum(x[i, j] for j in vertices if j > i) == 2, "Degree(%s)" % i)

    conshdlr = TSPconshdlr(x)

    model.includeConshdlr(conshdlr, "TSP", "TSP subtour eliminator",
                          chckpriority = -10, needscons = False)
    model.setBoolParam("misc/allowdualreds", False)

    model.setObjective(quicksum(distance[i, j] * x[i, j] for (i, j) in pairs(vertices)),
                       "minimize")

    return model, x

def gpp(V,E):
    """gpp -- model for the graph partitioning problem
    Parameters:
        - V: set/list of nodes in the graph
        - E: set/list of edges in the graph
    Returns a model, ready to be solved.
    """
    model = Model("gpp")

    x = {}
    y = {}
    for i in V:
        x[i] = model.addVar(vtype="B", name="x(%s)"%i)
    for (i,j) in E:
        y[i,j] = model.addVar(vtype="B", name="y(%s,%s)"%(i,j))

    model.addCons(quicksum(x[i] for i in V) == len(V)/2, "Partition")

    for (i,j) in E:
        model.addCons(x[i] - x[j] <= y[i,j], "Edge(%s,%s)"%(i,j))
        model.addCons(x[j] - x[i] <= y[i,j], "Edge(%s,%s)"%(j,i))

    model.setObjective(quicksum(y[i,j] for (i,j) in E), "minimize")

    model.data = x
    return model

def vrp(V, c, m, q, Q):
    """solve_vrp -- solve the vehicle routing problem.
       - start with assignment model (depot has a special status)
       - add cuts until all components of the graph are connected
    Parameters:
        - V: set/list of nodes in the graph
        - c[i,j]: cost for traversing edge (i,j)
        - m: number of vehicles available
        - q[i]: demand for customer i
        - Q: vehicle capacity
    Returns the optimum objective value and the list of edges used.
    """

    model = Model("vrp")
    vrp_conshdlr = VRPconshdlr()

    x = {}
    for i in V:
        for j in V:
            if j > i and i == V[0]:       # depot
                x[i,j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)"%(i,j))
            elif j > i:
                x[i,j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)"%(i,j))

    model.addCons(quicksum(x[V[0],j] for j in V[1:]) == 2*m, "DegreeDepot")
    for i in V[1:]:
        model.addCons(quicksum(x[j,i] for j in V if j < i) +
                        quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)

    model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j>i), "minimize")
    model.data = x

    return model, vrp_conshdlr

def test_instance(instance):
    s = Model()
    s.hideOutput()
    s.readProblem(instance)
    s.optimize()
    name = os.path.split(instance)[1]
    if name.rsplit('.',1)[1].lower() == 'gz':
        name = name.rsplit('.',2)[0]
    else:
        name = name.rsplit('.',1)[0]

    # we do not need the solution status
    primalbound = s.getObjVal()
    dualbound = s.getDualbound()

    # get solution data from solu file
    primalsolu = primalsolutions.get(name, None)
    dualsolu = dualsolutions.get(name, None)

    if s.getObjectiveSense() == 'minimize':
        assert relGE(primalbound, dualsolu)
        assert relLE(dualbound, primalsolu)
    else:
        if( primalsolu == infinity ): primalsolu = -infinity
        if( dualsolu == infinity ): dualsolu = -infinity
        assert relLE(primalbound, dualsolu)
        assert relGE(dualbound, primalsolu)

def gcp(V,E,K):
    """gcp -- model for minimizing the number of colors in a graph
    Parameters:
        - V: set/list of nodes in the graph
        - E: set/list of edges in the graph
        - K: upper bound on the number of colors
    Returns a model, ready to be solved.
    """
    model = Model("gcp")
    x,y = {},{}
    for k in range(K):
        y[k] = model.addVar(vtype="B", name="y(%s)"%k)
        for i in V:
            x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k))

    for i in V:
        model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)"%i)

    for (i,j) in E:
        for k in range(K):
            model.addCons(x[i,k] + x[j,k] <= y[k], "NotSameColor(%s,%s,%s)"%(i,j,k))

    model.setObjective(quicksum(y[k] for k in range(K)), "minimize")

    model.data = x
    return model

def mctransp(I,J,K,c,d,M):
    """mctransp -- model for solving the Multi-commodity Transportation Problem
    Parameters:
        - I: set of customers
        - J: set of facilities
        - K: set of commodities
        - c[i,j,k]: unit transportation cost on arc (i,j) for commodity k
        - d[i][k]: demand for commodity k at node i
        - M[j]: capacity
    Returns a model, ready to be solved.
    """

    model = Model("multi-commodity transportation")

    # Create variables
    x = {}

    for (i,j,k) in c:
        x[i,j,k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i,j,k))

    # Demand constraints
    for i in I:
        for k in K:
            model.addCons(sum(x[i,j,k] for j in J if (i,j,k) in x) == d[i,k], "Demand(%s,%s)" % (i,k))

    # Capacity constraints
    for j in J:
        model.addCons(sum(x[i,j,k] for (i,j2,k) in x if j2 == j) <= M[j], "Capacity(%s)" % j)

    # Objective
    model.setObjective(quicksum(c[i,j,k]*x[i,j,k]  for (i,j,k) in x), "minimize")

    model.data = x

    return model

def scheduling_time_index(J,p,r,w):
    """
    scheduling_time_index: model for the one machine total weighted tardiness problem

    Model for the one machine total weighted tardiness problem
    using the time index formulation

    Parameters:
        - J: set of jobs
        - p[j]: processing time of job j
        - r[j]: earliest start time of job j
        - w[j]: weighted of job j; the objective is the sum of the weighted completion time

    Returns a model, ready to be solved.
    """
    model = Model("scheduling: time index")
    T = max(r.values()) + sum(p.values())
    X = {}   # X[j,t]=1 if job j starts processing at time t, 0 otherwise
    for j in J:
        for t in range(r[j], T-p[j]+2):
            X[j,t] = model.addVar(vtype="B", name="x(%s,%s)"%(j,t))

    for j in J:
        model.addCons(quicksum(X[j,t] for t in range(1,T+1) if (j,t) in X) == 1, "JobExecution(%s)"%(j))

    for t in range(1,T+1):
        ind = [(j,t2) for j in J for t2 in range(t-p[j]+1,t+1) if (j,t2) in X]
        if ind != []:
            model.addCons(quicksum(X[j,t2] for (j,t2) in ind) <= 1, "MachineUB(%s)"%t)

    model.setObjective(quicksum((w[j] * (t - 1 + p[j])) * X[j,t] for (j,t) in X), "minimize")

    model.data = X
    return model

def eoq_soco(I,F,h,d,w,W):
    """eoq_soco --  multi-item capacitated economic ordering quantity model using soco
    Parameters:
        - I: set of items
        - F[i]: ordering cost for item i
        - h[i]: holding cost for item i
        - d[i]: demand for item i
        - w[i]: unit weight for item i
        - W: capacity (limit on order quantity)
    Returns a model, ready to be solved.
    """
    model = Model("EOQ model using SOCO")

    T,c = {},{}
    for i in I:
        T[i] = model.addVar(vtype="C", name="T(%s)"%i)  # cycle time for item i
        c[i] = model.addVar(vtype="C", name="c(%s)"%i)  # total cost for item i

    for i in I:
        model.addCons(F[i] <= c[i]*T[i])

    model.addCons(quicksum(w[i]*d[i]*T[i] for i in I) <= W)

    model.setObjective(quicksum(c[i] + h[i]*d[i]*T[i]*0.5 for i in I), "minimize")

    model.data = T,c
    return model

def test_quicksum_model():
    m = Model("quicksum")
    x = m.addVar("x")
    y = m.addVar("y")
    z = m.addVar("z")
    c = 2.3

    q = quicksum([x,y,z,c]) == 0.0
    s =      sum([x,y,z,c]) == 0.0

    assert(q.expr.terms == s.expr.terms)

def test_lpi():
    # create LP instance
    myLP = LP()

    myLP.addRow(entries = [(0,1),(1,2)] ,lhs = 5)
    lhs, rhs = myLP.getSides()
    assert lhs[0] == 5.0
    assert rhs[0] == myLP.infinity()

    assert(myLP.ncols() == 2)
    myLP.chgObj(0, 1.0)
    myLP.chgObj(1, 2.0)

    solval = myLP.solve()
    # create solver instance

    s = Model()
    # add some variables
    x = s.addVar("x", obj=1.0)
    y = s.addVar("y", obj=2.0)
    # add some constraint
    s.addCons(x + 2*y >= 5)
    # solve problem
    s.optimize()

    # check solution
    assert round(s.getVal(x)) == 5.0
    assert round(s.getVal(y)) == 0.0

    assert round(s.getObjVal() == solval)

def test_largequadratic():
    # inspired from performance issue on
    # http://stackoverflow.com/questions/38434300

    m = Model("dense_quadratic")
    dim = 200
    x = [m.addVar("x_%d" % i) for i in range(dim)]
    expr = quicksum((i+j+1)*x[i]*x[j]
                    for i in range(dim)
                    for j in range(dim))
    cons = expr <= 1.0
    #                              upper triangle,     diagonal
    assert len(cons.expr.terms) == dim * (dim-1) / 2 + dim
    m.addCons(cons)

def diet(F,N,a,b,c,d):
    """diet -- model for the modern diet problem
    Parameters:
        F - set of foods
        N - set of nutrients
        a[i] - minimum intake of nutrient i
        b[i] - maximum intake of nutrient i
        c[j] - cost of food j
        d[j][i] - amount of nutrient i in food j
    Returns a model, ready to be solved.
    """
    model = Model("modern diet")

    # Create variables
    x,y,z = {},{},{}
    for j in F:
        x[j] = model.addVar(vtype="I", name="x(%s)" % j)

    for i in N:
        z[i] = model.addVar(lb=a[i], ub=b[i], vtype="C", name="z(%s)" % i)

    # Constraints:
    for i in N:
        model.addCons(quicksum(d[j][i]*x[j] for j in F) == z[i], name="Nutr(%s)" % i)

    model.setObjective(quicksum(c[j]*x[j]  for j in F), "minimize")

    model.data = x,y,z
    return model

def flp_nonlinear_sos(I,J,d,M,f,c,K):
    """flp_nonlinear_sos --  use model with SOS constraints
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for customer i
        - M[j]: capacity of facility j
        - f[j]: fixed cost for using a facility in point j
        - c[i,j]: unit cost of servicing demand point i from facility j
        - K: number of linear pieces for approximation of non-linear cost function
    Returns a model, ready to be solved.
    """
    a,b = {},{}
    for j in J:
        U = M[j]
        L = 0
        width = U/float(K)
        a[j] = [k*width for k in range(K+1)]
        b[j] = [f[j]*math.sqrt(value) for value in a[j]]

    model = Model("nonlinear flp -- use model with SOS constraints")
    
    x = {}
    for j in J:
        for i in I:
            x[i,j] = model.addVar(vtype="C", name="x(%s,%s)"%(i,j))  # i's demand satisfied from j
    
    # total volume transported from plant j, corresponding (linearized) cost, selection variable:
    X,F,z = {},{},{}
    for j in J:
        # add constraints for linking piecewise linear part:
        X[j],F[j],z[j] = convex_comb_sos(model,a[j],b[j])
        X[j].ub = M[j]
        # for i in I:
        #     model.addCons(
        #         x[i,j] <= \
        #         quicksum(min(d[i],a[j][k+1]) * (z[j][k] + z[j][k+1])\
        #                  for k in range(len(a[j])-1)),
        #         "Strong(%s,%s)"%(i,j))

    # constraints for customer's demand satisfaction
    for i in I:
        model.addCons(quicksum(x[i,j] for j in J) == d[i], "Demand(%s)"%i)

    for j in J:
        model.addCons(quicksum(x[i,j] for i in I) == X[j], "Capacity(%s)"%j)

    model.setObjective(quicksum(F[j] for j in J) +\
                       quicksum(c[i,j]*x[i,j] for j in J for i in I),\
                       "minimize")
    
    model.data = x,X,F
    return model

def mctransp(I,J,K,c,d,M):
    """mctransp -- model for solving the Multi-commodity Transportation Problem
    Parameters:
        - I: set of customers
        - J: set of facilities
        - K: set of commodities
        - c[i,j,k]: unit transportation cost on arc (i,j) for commodity k
        - d[i][k]: demand for commodity k at node i
        - M[j]: capacity
    Returns a model, ready to be solved.
    """

    model = Model("multi-commodity transportation")

    # Create variables
    x = {}
    for (i,j,k) in c:
        x[i,j,k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i,j,k), obj=c[i,j,k])
# todo
    arcs = tuplelist([(i,j,k) for (i,j,k) in x])

    # Demand constraints
    for i in I:
        for k in K:
            model.addCons(sum(x[i,j,k] for (i,j,k) in arcs.select(i,"*",k)) == d[i,k], "Demand(%s,%s)" % (i,k))

    # Capacity constraints
    for j in J:
        model.addConstr(sum(x[i,j,k] for (i,j,k) in arcs.select("*",j,"*")) <= M[j], "Capacity(%s)" % j)

    model.data = x
    return model

def tsp(V,c):
    """tsp -- model for solving the traveling salesman problem with callbacks
       - start with assignment model
       - add cuts until there are no sub-cycles
    Parameters:
        - V: set/list of nodes in the graph
        - c[i,j]: cost for traversing edge (i,j)
    Returns the optimum objective value and the list of edges used.
    """
    model = Model("TSP_lazy")
    conshdlr = TSPconshdlr()

    x = {}
    for i in V:
        for j in V:
            if j > i:
                x[i,j] = model.addVar(vtype = "B",name = "x(%s,%s)" % (i,j))

    for i in V:
        model.addCons(quicksum(x[j, i] for j in V if j < i) +
                      quicksum(x[i, j] for j in V if j > i) == 2, "Degree(%s)" % i)

    model.setObjective(quicksum(c[i, j] * x[i, j] for i in V for j in V if j > i), "minimize")

    model.data = x
    return model, conshdlr

def test_lp():
    # create solver instance
    s = Model()

    # add some variables
    x = s.addVar("x", vtype='C', obj=1.0)
    y = s.addVar("y", vtype='C', obj=2.0)

    # add some constraint
    s.addCons(x + 2*y >= 5.0)

    # solve problem
    s.optimize()

    solution = s.getBestSol()

    # print solution
    assert (s.getVal(x) == s.getSolVal(solution, x))
    assert (s.getVal(y) == s.getSolVal(solution, y))
    assert round(s.getVal(x)) == 5.0
    assert round(s.getVal(y)) == 0.0

def test_quad_coeffs():
    """test coefficient access method for quadratic constraints"""
    scip = Model()
    x = scip.addVar()
    y = scip.addVar()
    z = scip.addVar()

    c = scip.addCons(2*x*y + 0.5*x**2 + 4*z >= 10)
    assert c.isQuadratic()

    bilinterms, quadterms, linterms = scip.getTermsQuadratic(c)

    assert bilinterms[0][0].name == x.name
    assert bilinterms[0][1].name == y.name
    assert bilinterms[0][2] == 2

    assert quadterms[0][0].name == x.name
    assert quadterms[0][1] == 0.5

    assert linterms[0][0].name == z.name
    assert linterms[0][1] == 4