# cogirth-packing gap of projections of graphic matroids
# see if gap(projection) <= 2 * gap(graph)
# this is wrong. see the last matrix in projection.out

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 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
f = lambda g: g.is_connected()
for N in range(4,10):
    for g in filter(f, graphs(N)):
        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
        V = VectorSpace(GF(2), N)
        for v in filter(lambda v: v.hamming_weight() % 2 == 0 and v!=0, V):
            cnt=cnt+1
            v_col=matrix(v).transpose()
            A_t=A.augment(v_col)
            # print(A_t)
            M=Matroid(matrix=A_t,field=GF(2))/m     #contract the last element
            # print(M)
            cogirth = base_hittingset_with_callback(M,integral=true)
            strength = base_hittingset_with_callback(M,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"n={n:<2}, m={m:<4}    graph gap = {fr(g_gap):8}   projection gap = {fr(m_gap):8}  factor={fr(factor):8}   max={maxfactor}")
        if factor >= 2:
            with open("projection.out", "a") as file:
                file.write(f"n={n:<2}, m={m:<4}    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)