solveCuttingStock: use column generation (Gilmore-Gomory approach).
    Parameters:
        - w: list of item's widths
        - q: number of items of a width
        - B: bin/roll capacity
    Returns a solution: list of lists, each of which with the cuts of a roll.

def solveCuttingStock(w,q,B):
    """solveCuttingStock: use column generation (Gilmore-Gomory approach).
    Parameters:
        - w: list of item's widths
        - q: number of items of a width
        - B: bin/roll capacity
    Returns a solution: list of lists, each of which with the cuts of a roll.
    """
    t = []      # patterns
    m = len(w)

    # Generate initial patterns with one size for each item width
    for (i,width) in enumerate(w):
        pat = [0]*m  # vector of number of orders to be packed into one roll (bin)
        pat[i] = int(B/width)
        t.append(pat)

    # if LOG:
    #     print "sizes of orders=",w
    #     print "quantities of orders=",q
    #     print "roll size=",B
    #     print "initial patterns",t

    K = len(t)
    master = Model("master LP") # master LP problem
    x = {}

    for k in range(K):
        x[k] = master.addVar(vtype="I", name="x(%s)"%k)

    orders = {}

    for i in range(m):
        orders[i] = master.addCons(
            quicksum(t[k][i]*x[k] for k in range(K) if t[k][i] > 0) >= q[i], "Order(%s)"%i)

    master.setObjective(quicksum(x[k] for k in range(K)), "minimize")

    # master.Params.OutputFlag = 0 # silent mode

    # iter = 0
    while True:
        # print "current patterns:"
        # for ti in t:
        #     print ti
        # print

        # iter += 1
        relax = master.relax()
        relax.optimize()
        pi = [relax.getDualsolLinear(c) for c in relax.getConss()] # keep dual variables

        knapsack = Model("KP")     # knapsack sub-problem
        knapsack.setMaximize       # maximize
        y = {}

        for i in range(m):
            y[i] = knapsack.addVar(lb=0, ub=q[i], vtype="I", name="y(%s)"%i)

        knapsack.addCons(quicksum(w[i]*y[i] for i in range(m)) <= B, "Width")

        knapsack.setObjective(quicksum(pi[i]*y[i] for i in range(m)), "maximize")

        knapsack.hideOutput() # silent mode
        knapsack.optimize()
        # if LOG:
        #     print "objective of knapsack problem:", knapsack.ObjVal
        if knapsack.getObjVal() < 1+EPS: # break if no more columns
            break

        pat = [int(y[i].X+0.5) for i in y]      # new pattern
        t.append(pat)
        # if LOG:
        #     print "shadow prices and new pattern:"
        #     for (i,d) in enumerate(pi):
        #         print "\t%5s%12s%7s" % (i,d,pat[i])
        #     print

        # add new column to the master problem
        col = Column()
        for i in range(m):
            if t[K][i] > 0:
                col.addTerms(t[K][i], orders[i])
        x[K] = master.addVar(obj=1, vtype="I", name="x(%s)"%K, column=col)

        # master.write("MP" + str(iter) + ".lp")
        K += 1


    # Finally, solve the IP
    # if LOG:
    #     master.Params.OutputFlag = 1 # verbose mode
    master.optimize()

    # if LOG:
    #     print
    #     print "final solution (integer master problem):  objective =", master.ObjVal
    #     print "patterns:"
    #     for k in x:
    #         if x[k].X > EPS:
    #             print "pattern",k,
    #             print "\tsizes:",
    #             print [w[i] for i in range(m) if t[k][i]>0 for j in range(t[k][i]) ],
    #             print "--> %s rolls" % int(x[k].X+.5)

    rolls = []
    for k in x:
        for j in range(int(x[k].X + .5)):
            rolls.append(sorted([w[i] for i in range(m) if t[k][i]>0 for j in range(t[k][i])]))
    rolls.sort()
    return rolls