def kcenter(I,J,c,k):
    """kcenter -- minimize the maximum travel cost from customers to k facilities.
    Parameters:
        - I: set of customers
        - J: set of potential facilities
        - c[i,j]: cost of servicing customer i from facility j
        - k: number of facilities to be used
    Returns a model, ready to be solved.
    """

    model = Model("k-center")
    z = model.addVar(vtype="C", name="z")
    x,y = {},{}

    for j in J:
        y[j] = model.addVar(vtype="B", name="y(%s)"%j)
        for i in I:
            x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))


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

        for j in J:
            model.addCons(x[i,j] <= y[j], "Strong(%s,%s)"%(i,j))
            model.addCons(c[i,j]*x[i,j] <= z, "Max_x(%s,%s)"%(i,j))

    model.addCons(quicksum(y[j] for j in J) == k, "Facilities")

    model.setObjective(z, "minimize")
    model.data = x,y

    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 p_portfolio(I,sigma,r,alpha,beta):
    """p_portfolio -- modified markowitz model for portfolio optimization.
    Parameters:
        - I: set of items
        - sigma[i]: standard deviation of item i
        - r[i]: revenue of item i
        - alpha: acceptance threshold
        - beta: desired confidence level
    Returns a model, ready to be solved.
    """

    model = Model("p_portfolio")

    x = {}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)  # quantity of i to buy
    rho = model.addVar(vtype="C", name="rho")
    rhoaux = model.addVar(vtype="C", name="rhoaux")

    model.addCons(rho == quicksum(r[i]*x[i] for i in I))
    model.addCons(quicksum(x[i] for i in I) == 1)

    model.addCons(rhoaux == (alpha - rho)*(1/phi_inv(beta))) #todo
    model.addCons(quicksum(sigma[i]**2 * x[i] * x[i] for i in I) <=  rhoaux * rhoaux)

    model.setObjective(rho, "maximize")

    model.data = x
    return model

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 prodmix(I,K,a,p,epsilon,LB):
    """prodmix:  robust production planning using soco
    Parameters:
        I - set of materials
        K - set of components
        a[i][k] -  coef. matrix
        p[i] - price of material i
        LB[k] - amount needed for k
    Returns a model, ready to be solved.
    """

    model = Model("robust product mix")

    x,rhs = {},{}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)
    for k in K:
        rhs[k] = model.addVar(vtype="C", name="rhs(%s)"%k)

    model.addCons(quicksum(x[i] for i in I) == 1)
    for k in K:
        model.addCons(rhs[k] == -LB[k]+ quicksum(a[i,k]*x[i] for i in I) )
        model.addCons(quicksum(epsilon*epsilon*x[i]*x[i] for i in I) <= rhs[k]*rhs[k])

    model.setObjective(quicksum(p[i]*x[i] for i in I), "minimize")

    model.data = x,rhs
    return model

def maxflow(V,M,source,sink):
    """maxflow: maximize flow from source to sink, taking into account arc capacities M
    Parameters:
        - V: set of vertices
        - M[i,j]: dictionary or capacity for arcs (i,j)
        - source: flow origin
        - sink: flow target
    Returns a model, ready to be solved.
        """
    # create max-flow underlying model, on which to find cuts
    model = Model("maxflow")
    
    f = {} # flow variable
    for (i,j) in M:
        f[i,j] = model.addVar(lb=-M[i,j], ub=M[i,j], name="flow(%s,%s)"%(i,j))

    cons = {}
    for i in V:
        if i != source and i != sink:
            cons[i] = model.addCons(
                quicksum(f[i,j] for j in V if i<j and (i,j) in M) - \
                quicksum(f[j,i] for j in V if i>j and (j,i) in M) == 0,
                "FlowCons(%s)"%i)

    model.setObjective(quicksum(f[i,j] for (i,j) in M if i==source), "maximize")

    # model.write("tmp.lp")
    model.data = f,cons
    return model

def markowitz(I,sigma,r,alpha):
    """markowitz -- simple markowitz model for portfolio optimization.
    Parameters:
        - I: set of items
        - sigma[i]: standard deviation of item i
        - r[i]: revenue of item i
        - alpha: acceptance threshold
    Returns a model, ready to be solved.
    """
    model = Model("markowitz")

    x = {}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)  # quantity of i to buy

    model.addCons(quicksum(r[i]*x[i] for i in I) >= alpha)
    model.addCons(quicksum(x[i] for i in I) == 1)

    # set nonlinear objective: SCIP only allow for linear objectives hence the following
    obj = model.addVar(vtype="C", name="objective", lb = None, ub = None)  # auxiliary variable to represent objective
    model.addCons(quicksum(sigma[i]**2 * x[i] * x[i] for i in I) <= obj)
    model.setObjective(obj, "minimize")

    model.data = x
    return model

    def addcut(self, checkonly, sol):
        D,Ts = self.data
        y,x,I = self.model.data
        cutsadded = False

        for ell in Ts:
            lhs = 0
            S,L = [],[]
            for t in range(1,ell+1):
                yt = self.model.getSolVal(sol, y[t])
                xt = self.model.getSolVal(sol, x[t])
                if D[t,ell]*yt < xt:
                    S.append(t)
                    lhs += D[t,ell]*yt
                else:
                    L.append(t)
                    lhs += xt
            if lhs < D[1,ell]:
                if checkonly:
                    return True
                else:
                    # add cutting plane constraint
                    self.model.addCons(quicksum([x[t] for t in L]) + \
                                       quicksum(D[t, ell] * y[t] for t in S)
                                       >= D[1, ell], removable = True)
                cutsadded = True
        return cutsadded

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 convex_comb_agg(model,a,b):
    """convex_comb_agg -- add piecewise relation convex combination formulation -- non-disaggregated.
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.
    """
    K = len(a)-1
    w,z = {},{}
    for k in range(K+1):
        w[k] = model.addVar(lb=0, ub=1, vtype="C")
    for k in range(K):
        z[k] = model.addVar(vtype="B")
    X = model.addVar(lb=a[0], ub=a[K], vtype="C")
    Y = model.addVar(lb=-model.infinity(), vtype="C")

    model.addCons(X == quicksum(a[k]*w[k] for k in range(K+1)))
    model.addCons(Y == quicksum(b[k]*w[k] for k in range(K+1)))
    model.addCons(w[0] <= z[0])
    model.addCons(w[K] <= z[K-1])
    for k in range(1,K):
        model.addCons(w[k] <= z[k-1]+z[k])
    model.addCons(quicksum(w[k] for k in range(K+1)) == 1)
    model.addCons(quicksum(z[k] for k in range(K)) == 1)
    return X,Y,z

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 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 mult_selection(model,a,b):
    """mult_selection -- add piecewise relation with multiple selection formulation
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.
    """

    K = len(a)-1
    w,z = {},{}
    for k in range(K):
        w[k] = model.addVar(lb=-model.infinity()) # do not name variables for avoiding clash
        z[k] = model.addVar(vtype="B")
    X = model.addVar(lb=a[0], ub=a[K], vtype="C")
    Y = model.addVar(lb=-model.infinity())

    for k in range(K):
        model.addCons(w[k] >= a[k]*z[k])
        model.addCons(w[k] <= a[k+1]*z[k])

    model.addCons(quicksum(z[k] for k in range(K)) == 1)
    model.addCons(X == quicksum(w[k] for k in range(K)))

    c = [float(b[k+1]-b[k])/(a[k+1]-a[k]) for k in range(K)]
    d = [b[k]-c[k]*a[k] for k in range(K)]
    model.addCons(Y == quicksum(d[k]*z[k] + c[k]*w[k] for k in range(K)))

    return X,Y,z

def convex_comb_dis(model,a,b):
    """convex_comb_dis -- add piecewise relation with convex combination formulation
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.
    """
    K = len(a)-1
    wL,wR,z = {},{},{}
    for k in range(K):
        wL[k] = model.addVar(lb=0, ub=1, vtype="C")
        wR[k] = model.addVar(lb=0, ub=1, vtype="C")
        z[k] = model.addVar(vtype="B")
    X = model.addVar(lb=a[0], ub=a[K], vtype="C")
    Y = model.addVar(lb=-model.infinity(), vtype="C")

    model.addCons(X == quicksum(a[k]*wL[k] + a[k+1]*wR[k] for k in range(K)))
    model.addCons(Y == quicksum(b[k]*wL[k] + b[k+1]*wR[k] for k in range(K)))
    for k in range(K):
        model.addCons(wL[k] + wR[k] == z[k])

    model.addCons(quicksum(z[k] for k in range(K)) == 1)

    return X,Y,z

def mils_echelon(T,K,P,f,g,c,d,h,a,M,UB,phi):
    """
    mils_echelon: echelon formulation for the multi-item, multi-stage lot-sizing problem

    Parameters:
        - T: number of periods
        - K: set of resources
        - P: set of items
        - f[t,p]: set-up costs (on period t, for product p)
        - g[t,p]: set-up times
        - c[t,p]: variable costs
        - d[t,p]: demand values
        - h[t,p]: holding costs
        - a[t,k,p]: amount of resource k for producing p in period t
        - M[t,k]: resource k upper bound on period t
        - UB[t,p]: upper bound of production time of product p in period t
        - phi[(i,j)]: units of i required to produce a unit of j (j parent of i)
    """
    rho = calc_rho(phi) # rho[(i,j)]: units of i required to produce a unit of j (j ancestor of i)

    model = Model("multi-stage lotsizing -- echelon formulation")

    y,x,E,H = {},{},{},{}
    Ts = range(1,T+1)
    for p in P:
        for t in Ts:
            y[t,p] = model.addVar(vtype="B", name="y(%s,%s)"%(t,p))
            x[t,p] = model.addVar(vtype="C", name="x(%s,%s)"%(t,p))
            H[t,p] = h[t,p] - sum([h[t,q]*phi[q,p] for (q,p2) in phi if p2 == p])
            E[t,p] = model.addVar(vtype="C", name="E(%s,%s)"%(t,p))        # echelon inventory
        E[0,p] = model.addVar(vtype="C", name="E(%s,%s)"%(0,p))    # echelon inventory

    for t in Ts:
        for p in P:
            # flow conservation constraints
            dsum = d[t,p] + sum([rho[p,q]*d[t,q] for (p2,q) in rho if p2 == p])
            model.addCons(E[t-1,p] + x[t,p] == E[t,p] + dsum, "FlowCons(%s,%s)"%(t,p))

            # capacity connection constraints
            model.addCons(x[t,p] <= UB[t,p]*y[t,p], "ConstrUB(%s,%s)"%(t,p))

        # time capacity constraints
        for k in K:
            model.addCons(quicksum(a[t,k,p]*x[t,p] + g[t,p]*y[t,p] for p in P) <= M[t,k],
                            "TimeUB(%s,%s)"%(t,k))


    # calculate echelon quantities
    for p in P:
        model.addCons(E[0,p] == 0, "EchelonInit(%s)"%(p))
        for t in Ts:
            model.addCons(E[t,p] >= quicksum(phi[p,q]*E[t,q] for (p2,q) in phi if p2 == p),
                            "EchelonLB(%s,%s)"%(t,p))

    model.setObjective(\
        quicksum(f[t,p]*y[t,p] + c[t,p]*x[t,p] + H[t,p]*E[t,p] for t in Ts for p in P), \
        "minimize")

    model.data = y,x,E
    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 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 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 eld_another(U,p_min,p_max,d,brk):
    """eld -- economic load dispatching in electricity generation
    Parameters:
        - U: set of generators (units)
        - p_min[u]: minimum operating power for unit u
        - p_max[u]: maximum operating power for unit u
        - d: demand
        - brk[u][k]: (x,y) coordinates of breakpoint k, k=0,...,K for unit u
    Returns a model, ready to be solved.
    """
    model = Model("Economic load dispatching")

    # set objective based on piecewise linear approximation
    p,F,z = {},{},{}
    for u in U:
        abrk = [X for (X,Y) in brk[u]]
        bbrk = [Y for (X,Y) in brk[u]]
        p[u],F[u],z[u] = convex_comb_sos(model,abrk,bbrk)
        p[u].lb = p_min[u]
        p[u].ub = p_max[u]

    # demand satisfaction
    model.addCons(quicksum(p[u] for u in U) == d, "demand")

    # objective
    model.setObjective(quicksum(F[u] for u in U), "minimize")

    model.data = p
    return model

def mtz_strong(n,c):
    """mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
    (potential formulation, adding stronger constraints)
    Parameters:
        n - number of nodes
        c[i,j] - cost for traversing arc (i,j)
    Returns a model, ready to be solved.
    """

    model = Model("atsp - mtz-strong")
    
    x,u = {},{}
    for i in range(1,n+1):
        u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i)
        for j in range(1,n+1):
            if i != j:
                x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))

    for i in range(1,n+1):
        model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
        model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)

    for i in range(1,n+1):
        for j in range(2,n+1):
            if i != j:
                model.addCons(u[i] - u[j] + (n-1)*x[i,j] + (n-3)*x[j,i] <= n-2, "LiftedMTZ(%s,%s)"%(i,j))

    for i in range(2,n+1):
        model.addCons(-x[1,i] - u[i] + (n-3)*x[i,1] <= -2, name="LiftedLB(%s)"%i)
        model.addCons(-x[i,1] + u[i] + (n-3)*x[1,i] <= n-2, name="LiftedUB(%s)"%i)

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