# 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)