gap is not constant

.gitignore

@@ -1,3 +1,6 @@
+
+.vscode/
+*.sage.py
 ## Core latex/pdflatex auxiliary files:
 *.aux
 *.lof

test.sage

@@ -0,0 +1,99 @@
+'''
+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()