102 lines
3.8 KiB
Python
102 lines
3.8 KiB
Python
# 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) |