'''
gap is not constant.

OPT(IP): 3.0
OPT(LP): 0.6000000000000003
gap: 4.999999999999997
w: {(0, 2): 7, (0, 3): 6, (1, 2): 3, (1, 3): 9}
c: {(0, 2): 8, (0, 3): 9, (1, 2): 10, (1, 3): 3}
b: 8.0
edges: [(0, 2), (0, 3), (1, 2), (1, 3)]
IP solution: [('x[0,2]', 0.0), ('x[0,3]', 0.0), ('x[1,2]', 1.0), ('x[1,3]', 0.0), ('y[0,2]', 0.0), ('y[0,3]', 0.0), ('y[1,2]', 0.0), ('y[1,3]', 1.0)]
LP solution: [('x[0,2]', 0.0), ('x[0,3]', 0.0), ('x[1,2]', 0.20000000000000012), ('x[1,3]', 0.0), ('y[0,2]', 0.39999999999999997), ('y[0,3]', 0.39999999999999997), ('y[1,2]', 0.0), ('y[1,3]', 0.39999999999999997)]
'''

from random import randint

from sage.all import *
from sage.matroids.all import *

import gurobipy as gp
from gurobipy import GRB

def interdiction(G, w, c, b, vartype=GRB.BINARY):
    """
    Interdiction problem on a graph G with cut weights w, interdiction costs c 
    and budget b.
    """
    model = gp.Model("mip1")
    # model.Params.LogToConsole = 0
    edges=G.edges(labels=False)
    # Variables
    x = model.addVars(edges, vtype=vartype, name="x", lb=0)
    y = model.addVars(edges, vtype=vartype, name="y", lb=0)

    # Objective
    model.setObjective(gp.quicksum(w[e] * x[e] for e in edges), 
                       GRB.MINIMIZE)

    # Constraints
    # knapsack
    model.addConstr(gp.quicksum(c[e] * y[e] for e in edges) <= b)
    # cut
    M=Matroid(G)
    for b in M.bases():
        model.addConstr(gp.quicksum(x[e]+y[e] for e in b) >= 1)

    # Solve
    model.optimize()
    solutions = []
    for v in model.getVars():
        solutions.append((v.varName, v.x))

    return model.ObjVal,solutions

def mincut(G, c, vartype=GRB.BINARY):
    model = gp.Model("mip1")
    # model.Params.LogToConsole = 0
    edges=G.edges(labels=False)
    # Variables
    x = model.addVars(edges, vtype=vartype, name="x", lb=0)
    # Objective
    model.setObjective(gp.quicksum(c[e] * x[e] for e in edges), 
                       GRB.MINIMIZE)
    # Constraints
    M=Matroid(G)
    for b in M.bases():
        model.addConstr(gp.quicksum(x[e] for e in b) >= 1)
    # Solve
    model.optimize()
    return model.ObjVal

## enumerate all graphs
lb = 4
ub = 6

f = lambda G: G.is_connected()
for n in range(lb, ub + 1):
    for G in filter(f, graphs(n)):
        for _ in range(10): # use 10 random weights
            # generate w,c,b
            edges=G.edges(labels=False)
            w = {e: randint(1, 10) for e in edges}
            c = {e: randint(1, 10) for e in edges}
            b = mincut(G, c)-randint(1, 5)
            if b < 0: continue

            IP,solIP=interdiction(G, w, c, b, GRB.BINARY)
            LP,solLP=interdiction(G, w, c, b, GRB.CONTINUOUS)
            print(f"OPT(IP): {IP}")
            print(f"OPT(LP): {LP}")
            print(f"gap: {IP/LP}")
            if IP/LP > 4:
                print(f"w: {w}")
                print(f"c: {c}")
                print(f"b: {b}")
                print("edges:", G.edges(labels=False))
                print("IP solution:", solIP)
                print("LP solution:", solLP)
                exit()