testing my 2 factor conjecture

Kmn_2factor.sage

@@ -0,0 +1,102 @@
+# cogirth-packing gap of projections of graphic matroids
+# see if gap(projection) <= 2 * gap(graph)
+
+from sage.all import *
+from sage.matroids.all import *
+from sage.graphs.all import *
+import gurobipy as gp
+from gurobipy import GRB
+from fractions import Fraction
+
+env = gp.Env(empty=True)
+env.setParam("OutputFlag",0)
+env.start()
+
+def representative_vectors(m, n):
+    for w1 in range(m+1):
+        for w2 in range(n+1):
+            v = [0]*(m+n)
+            v[:w1] = [1]*w1
+            v[m:m+w2] = [1]*w2
+            yield tuple(v)
+
+def base_hittingset_with_callback(M, integral=true):
+    # model
+    model = gp.Model("mip1",env=env)
+    groundset = M.groundset()
+    x = dict()
+    if integral:
+        for e in groundset: x[e]=model.addVar(vtype=GRB.BINARY)
+    else:
+        for e in groundset: x[e]=model.addVar(vtype=GRB.CONTINUOUS,lb=0,ub=1)
+    # there is no lazy constraint for LP...
+    for B in M.bases():
+        model.addConstr(gp.quicksum([x[e] for e in B])>=1)
+    model.setObjective(gp.quicksum([x[e] for e in groundset]), GRB.MINIMIZE)
+    # callback
+    def callback_func(model, where):
+        if integral and where == GRB.Callback.MIPSOL:
+            sol_values = {key: model.cbGetSolution(var) 
+                          for key, var in x.items()}
+            # find min weight base in the matroid
+            base = frozenset()
+            minweight=0
+            for e in sorted(groundset, key=lambda ee: sol_values[ee]):
+                if M.is_independent(base.union([e])):
+                    base=base.union([e])
+                    minweight=minweight+sol_values[e]
+            if minweight < 1:
+                model.cbLazy(gp.quicksum([x[e] for e in base])>=1)
+    # solve
+    if integral:
+        model.Params.Heuristics = 0
+        model.Params.LazyConstraints = 1
+        model.optimize(callback_func)
+    else:
+        model.optimize()
+    return model.ObjVal
+
+
+cnt=0   # actual number of instances tested
+maxfactor=0
+for N in range(2,10):
+    for M in range(N,10):
+        g=graphs.CompleteBipartiteGraph(N,M)
+        MG=Matroid(g)
+        mincut=base_hittingset_with_callback(MG)
+        treepacking=base_hittingset_with_callback(MG,integral=false)
+        g_gap=mincut/treepacking
+        A=g.incidence_matrix()
+        n,m = A.dimensions()
+        # enumerate all vectors in F_2^n with even number of 1s
+        m_gap=0
+        for v in representative_vectors(N,M):
+            cnt=cnt+1
+            v_col=matrix(v).transpose()
+            A_t=A.augment(v_col)
+            # print(A_t)
+            MM=Matroid(matrix=A_t,field=GF(2))/m     #contract the last element
+            cogirth = base_hittingset_with_callback(MM,integral=true)
+            strength = base_hittingset_with_callback(MM,integral=false)
+            gap = cogirth/strength
+            m_gap=max(m_gap,gap)
+            # maxgap=max(gap,maxgap)
+            # if gap > maxgap:
+            #     maxgap = gap
+            #     print(f"find a large gap: {gap}")
+            #     with open("projection.out", "a") as file:
+            #         file.write("##################################\n"
+            #                    +str(gap)+"\n"+str(A_t)
+            #                    +"\n##################################\n")
+        factor=m_gap/g_gap
+        fr = lambda xx:  str(Fraction(xx).limit_denominator(m))
+        maxfactor=max(maxfactor,factor)
+        print(f"K_{{{N},{M}}}    graph gap = {fr(g_gap):8}   projection gap = {fr(m_gap):8}  factor={fr(factor):8}   max={maxfactor}")
+        if factor >= 2:
+            with open("Kmn_2factor.out", "a") as file:
+                file.write(f"K_{{{N},{M}}}    graph gap = {fr(g_gap):8}   projection gap = {fr(m_gap):8}  factor={fr(factor):8}   max={maxfactor}\n"
+                            +str(A_t)+"\n")
+        if factor > 2:
+            print("conjecture is wrong")
+            print(str(A_t))
+            exit(1)