"""
lotsizing.py:  solve the multi-item lot-sizing problem.

Approaches:
    - mils: solve the problem using the standard formulation
    - mils_fl: solve the problem using the facility location (tighten) formulation

Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
import random
from gurobipy import *

def mils(T,P,f,g,c,d,h,M):
    """
    mils: standard formulation for the multi-item lot-sizing problem
    Parameters:
        - T: number of periods
        - P: set of products
        - 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
        - M[t]: resource upper bound on period t
    Returns a model, ready to be solved.
    """
    def mils_callback(model,where):
        # remember to set     model.params.DualReductions = 0     before using!
        if where != GRB.Callback.MIPSOL and where != GRB.Callback.MIPNODE:
            return
        for p in P:
            for ell in Ts:
                lhs = 0
                S,L = [],[]
                for t in range(1,ell+1):
                    yt = model.cbGetSolution(y[t,p])
                    xt = model.cbGetSolution(x[t,p])
                    if D[t,ell,p]*yt < xt:
                        S.append(t)
                        lhs += D[t,ell,p]*yt
                    else:
                        L.append(t)
                        lhs += xt
                if lhs < D[1,ell,p]:
                    # add cutting plane constraint
                    model.cbLazy(quicksum(x[t,p] for t in L) +\
                                 quicksum(D[t,ell,p] * y[t,p] for t in S)
                                 >= D[1,ell,p])
        return

    model = Model("standard multi-item lotsizing")
    y,x,I = {},{},{}
    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))
            I[t,p] = model.addVar(vtype="C", name="I(%s,%s)"%(t,p))
        I[0,p] = 0
    model.update()

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

        for p in P:
            # flow conservation constraints
            model.addConstr(I[t-1,p] + x[t,p] == I[t,p] + d[t,p], "FlowCons(%s,%s)"%(t,p))

            # capacity connection constraints
            model.addConstr(x[t,p] <= (M[t]-g[t,p])*y[t,p], "ConstrUB(%s,%s)"%(t,p))

            # tighten constraints
            model.addConstr(x[t,p] <= d[t,p]*y[t,p] + I[t,p], "Tighten(%s,%s)"%(t,p))

    model.setObjective(\
        quicksum(f[t,p]*y[t,p] + c[t,p]*x[t,p] + h[t,p]*I[t,p] for t in Ts for p in P),\
        GRB.MINIMIZE)

    # compute D[i,j,p] = sum_{t=i}^j d[t,p]
    D = {}
    for p in P:
        for t in Ts:
            s = 0
            for j in range(t,T+1):
                s += d[j,p]
                D[t,j,p] = s

    model.__data = y,x,I
    return model,mils_callback


def mils_fl(T,P,f,g,c,d,h,M):
    """
    mils_fl: facility location formulation for the multi-item lot-sizing problem

    Requires more variables, but gives a better solution because LB is
    better than the standard formulation.  It can be used as a
    heuristic method that is sometimes better than relax-and-fix.

    Parameters:
        - T: number of periods
        - P: set of products
        - 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
        - M[t]:   resource upper bound on period t
    Returns a model, ready to be solved.
    """
    Ts = range(1,T+1)

    model = Model("multi-item lotsizing -- facility location formulation")

    y,X = {},{}
    for p in P:
        for t in Ts:
            y[t,p] = model.addVar(vtype="B", name="y(%s,%s)"%(t,p))
            for s in range(1,t+1):
                X[s,t,p] = model.addVar(name="X(%s,%s,%s)"%(s,t,p))
    model.update()

    for t in Ts:
        # capacity constraints
        model.addConstr(quicksum(X[t,s,p] for s in range(t,T+1) for p in P) + \
                        quicksum(g[t,p]*y[t,p] for p in P) <= M[t],
                        "Capacity(%s)"%(t))

        for p in P:
            # demand satisfaction constraints
            model.addConstr(quicksum(X[s,t,p] for s in range(1,t+1)) == d[t,p], "Demand(%s,%s)"%(t,p))
            
            # connection constraints
            for s in range(1,t+1):
                model.addConstr(X[s,t,p] <= d[t,p] * y[s,p], "Connect(%s,%s,%s)"%(s,t,p))
            
    C = {} # variable costs plus holding costs
    for p in P:
        for s in Ts:
            sumC = 0
            for t in range(s,T+1):
                C[s,t,p] = (c[s,p] + sumC)
                sumC += h[t,p]

    model.setObjective(quicksum(f[t,p]*y[t,p] for t in Ts for p in P) + \
                       quicksum(C[s,t,p]*X[s,t,p] for t in Ts for p in P for s in range(1,t+1)),
                       GRB.MINIMIZE)

    model.update()
    model.__data = y,X
    model.write("tmp.lp")
    return model



def trigeiro(T,N,factor):
    """
    Data generator for the multi-item lot-sizing problem
    it uses a simular algorithm for generating the standard benchmarks in:
    "Capacitated Lot Sizing with Setup Times" by
    William W. Trigeiro, L. Joseph Thomas, John O. McClain
    MANAGEMENT SCIENCE
    Vol. 35, No. 3, March 1989, pp. 353-366

    Parameters:
        - T: number of periods
        - N: number of products
        - factor: value for controlling constraining factor of capacity:
            - 0.75:  lightly-constrained instances
            - 1.10:  constrained instances
    """
    P = range(1,N+1)
    f,g,c,d,h,M = {},{},{},{},{},{}

    sumT = 0
    for t in range(1,T+1):
        for p in P:
            # capacity used per unit production: 1, except for 
            # except for specific instances with random value in [0.5, 1.5]
            # (not tackled in our model)

            # setup times
            g[t,p] = 10 * random.randint(1,5)   # 10, 50: trigeiro's values

            # set-up costs
            f[t,p] = 100 * random.randint(1,10) # checked from Wolsey's instances
            c[t,p] = 0                          # variable costs

            # demands
            d[t,p] = 100+random.randint(-25,25) # checked from Wolsey's instances
            if t <= 4:
                if random.random() < .25:       # trigeiro's parameter
                    d[t,p] = 0
            sumT += g[t,p] + d[t,p]             # sumT is the total capacity usage in the lot-for-lot solution
            h[t,p] = random.randint(1,5)        # holding costs; checked from Wolsey's instances

    for t in range(1,T+1):
        M[t] = int(float(sumT)/float(T)/factor)

    return P,f,g,c,d,h,M
    


if __name__ == "__main__":

    # test only a subset of instances
    T,N,factor = 15,6,0.75
    P,f,g,c,d,h,M = trigeiro(T,N,factor)
    print "\n\n\nStandard formulation + cutting plane ======================"
    model,mils_callback = mils(T,P,f,g,c,d,h,M)
    model.params.DualReductions = 0
    model.optimize(mils_callback)
    y,x,I = model.__data
 
    print "Opt.value=",model.ObjVal
    standard = model.ObjVal
    print "%7s%7s%7s%7s%12s%12s" % ("prod","t","dem","y","x","I")
    for p in P:
        for t in range(1,T+1):
            print "%7s%7s%7s%7s%12s%12s" % \
                  (p,t,d[t,p],int(y[t,p].X+0.5),round(x[t,p].X,5),round(I[t,p].X,5))

    exit(0)
    # to test all instances: remove 'exit' above

    # sizes: 4 to 36 items, 15 to 30 time periods;   table 2: T = 15, 30, N = 6, 12, 24
    # all costs constant over time
    for T in [15,30]:
        for N in [6,12,24]:
            for factor in [0.75,1.0,1.1]:
                # 0.75: lightly-constrained instances
                # 1.1:  constrained instances

                random.seed(1)
                # T,P = 10,5   # number of periods and number of products
                # # T,P = 30,24  # number of periods and number of products
                # P,f,g,c,d,h,M = make_data(T,P)
                P,f,g,c,d,h,M = trigeiro(T,N,factor)
                #for i in d:
                #    d[i] = 100
                # print "f",f
                # print "g",g
                # print "c",c
                # print "d",d
                # print "h",h
                # print "M",M
                print "\n\n\n*** Trigeiros instance for %s periods, %s products, %s capacity factor" %\
                      (T,N,factor)

                # standard formulation
                print "\n\n\nStandard formulation======================"
                model,_ = mils(T,P,f,g,c,d,h,M)
                model.optimize()
                y,x,I = model.__data
             
                status = model.Status
                if status == GRB.Status.UNBOUNDED:
                    print "Unbounded instance"
                elif status == GRB.Status.OPTIMAL:
                    print "Opt.value=",model.ObjVal
                    standard = model.ObjVal
                    print "%7s%7s%7s%7s%12s%12s" % ("prod","t","dem","y","x","I")
                    for p in P:
                        for t in range(1,T+1):
                            print "%7s%7s%7s%7s%12s%12s" % \
                                  (p,t,d[t,p],int(y[t,p].X+0.5),round(x[t,p].X,5),round(I[t,p].X,5))
                elif status != GRB.Status.INF_OR_UNBD and status != GRB.Status.INFEASIBLE:
                    print "Optimization stopped with status",status
                else:
                    print "Instance infeasible; computing IIS"
                    model.computeIIS()
                    print "\nThe following constraint(s) cannot be satisfied:"
                    for i in model.getConstrs():
                      if i.IISConstr:
                        print i.ConstrName
             
             
                # standard formulation + cutting plane
                print "\n\n\nStandard formulation + cutting plane ======================"
                model,mils_callback = mils(T,P,f,g,c,d,h,M)
                model.params.DualReductions = 0
                model.optimize(mils_callback)
                y,x,I = model.__data
             
                status = model.Status
                if status == GRB.Status.UNBOUNDED:
                    print "Unbounded instance"
                elif status == GRB.Status.OPTIMAL:
                    print "Opt.value=",model.ObjVal
                    standard = model.ObjVal
                    print "%7s%7s%7s%7s%12s%12s" % ("prod","t","dem","y","x","I")
                    for p in P:
                        for t in range(1,T+1):
                            print "%7s%7s%7s%7s%12s%12s" % \
                                  (p,t,d[t,p],int(y[t,p].X+0.5),round(x[t,p].X,5),round(I[t,p].X,5))
                    assert abs(model.objval - standard) < 1.e-6
                elif status != GRB.Status.INF_OR_UNBD and status != GRB.Status.INFEASIBLE:
                    print "Optimization stopped with status",status
                else:
                    print "Instance infeasible; computing IIS"
                    model.computeIIS()
                    print "\nThe following constraint(s) cannot be satisfied:"
                    for i in model.getConstrs():
                      if i.IISConstr:
                        print i.ConstrName



                #  facility location formulation
                print "\n\n\nFacility location formulation======================"
                model = mils_fl(T,P,f,g,c,d,h,M)
                model.optimize()
                y,X = model.__data
             
                status = model.Status
                if status == GRB.Status.UNBOUNDED:
                    print "Unbounded instance"
                elif status == GRB.Status.OPTIMAL:
                    print "Opt.value=",model.ObjVal
                    print "%7s%7s%7s%7s%12s%12s" % ("prod","t","dem","y","sum(X)","inv.mov")
                    for p in P:
                        for t in range(1,T+1):
                            xd = sum(X[t,s,p].X for s in range(t,T+1))
                            I = xd - d[t,p]
                            print "%7s%7s%7s%7s%12s%12s" % \
                                  (p,t,d[t,p],int(y[t,p].X+0.5),round(xd,5),round(I))
                    assert abs(model.objval - standard) < 1.e-6
                elif status != GRB.Status.INF_OR_UNBD and status != GRB.Status.INFEASIBLE:
                    print "Optimization stopped with status",status
                else:
                    print "Instance infeasible; computing IIS"
                    model.computeIIS()
                    print "\nThe following constraint(s) cannot be satisfied:"
                    for i in model.getConstrs():
                      if i.IISConstr:
                        print i.ConstrName