# Cutting stock model with column generation.
#
# (C) 2020-2025 Fair Isaac Corporation

import xpress as xp
import numpy as np

# Problem data:
# widths = [17.0, 21.0, 22.5, 24.0, 29.5]
# demand = [150, 96, 48, 108, 227]
# max_width = 94.0
widths = [6, 11, 17, 21, 24, 28, 30, 33, 42, 49, 56, 69, 74, 87, 91]
demand = [9, 6, 20, 30, 17, 19, 25, 12, 8, 20, 5, 14, 15, 18, 10]
max_width = 100
eps = 1e-6
max_cols = 100

def solve_knapsack_problem(n_item_types, profit, resource_consumed, total_resource, demand):
    """
    Solve knapsack subproblem.
    :param n_item_types: Total number of items
    :param profit: Profit of each item. 1-d numpy array
    :param resource_consumed: Resource consumed by each item. 1-d numpy array
    :param total_resource: total available resource
    :param demand: Demand for each item. 1-d numpy array
    :return:
        xbest: Number of items of type i used in an optimal solution
        zbest: Value of the optimal solution
    """
    N = range(n_item_types)
    p = xp.problem(name="knapsack")
    p.controls.outputlog = 0

    # Create a NumPy array of variables
    x = xp.array([p.addVariable(name='use_{}'.format(i), lb=0, ub=demand[i], vartype=xp.integer) for i in N])

    p.addConstraint(xp.Dot(resource_consumed, x) <= total_resource)
    p.setObjective(xp.Dot(profit, x), sense=xp.maximize)
    p.optimize()
    return p.attributes.objval, p.getSolution()

def cutting_stock_model():
    """
    Solve a cutting stock problem by defining main problem and
    subproblem.
    """
    n_widths = len(widths)
    N = range(n_widths)
    p = xp.problem(name='CutStock')
    patterns = np.zeros((n_widths, n_widths))
    for j in N:
        patterns[j][j] = int(max_width / widths[j])

    # Create array of variables
    v_patterns = xp.array([
        p.addVariable(name='pat_{}'.format(i), lb=0, ub=int(float(demand[i] / patterns[i][i]) + 1))
        for i in N
    ])
    p.setObjective(xp.Dot(np.ones(n_widths), v_patterns), sense=xp.minimize)
    p.addConstraint(xp.Dot(patterns, v_patterns) >= demand)
    p.controls.outputlog = 0

    for npass in range(max_cols):
        print('=============Iteration {}==========='.format(npass))
        p.optimize()
        cs_lp_sol = p.getSolution()
        duals = p.getDuals()
        obj_val = p.attributes.objval
        z, x_best = solve_knapsack_problem(n_widths, np.array(duals),
            np.array(widths), max_width, demand)
        if z < 1 + eps:
            print('No profitable column found')
            break
        print('New pattern found with marginal cost {}'.format(z - 1));
        print(x_best)
        cur_pat = p.addVariable(name='pat_{}'.format(n_widths + npass),
            vartype=xp.continuous)

        for i in N:
            p.chgCoef(i, cur_pat, x_best[i])
        p.chgObj([cur_pat], [1.0])

    print("LP Solution after column generation. Obj: {}, Sol = {}".format(obj_val, cs_lp_sol))
    n_final_patterns = p.getAttrib('cols')
    print('Number of patterns in the final LP = {}'.format(n_final_patterns))
    # Change the variables to integers.
    p.chgColType([i for i in range(n_final_patterns)], ['I' for _ in range(n_final_patterns)])
    p.optimize()
    integer_obj_val = p.attributes.objval
    int_sol = p.getSolution()
    print("MIP Solution with final columns. Obj: {}, Sol = {}".format(integer_obj_val, int_sol))

if __name__ == '__main__':
    cutting_stock_model()