#### an optimized version of clique ####

def S(l):
    if (l[0] < l[1]):
        return (l[0],l[1]);
    else:
        return (l[1],l[0]);

def clique_minor(G, H, solver=None, verbose=0):
    
    from sage.numerical.mip import MixedIntegerLinearProgram, MIPSolverException, Sum
    from sage.rings.integer import Integer;
    p = MixedIntegerLinearProgram(solver=solver)
    p.solver_parameter("CPX_PARAM_THREADS", 1)
    
    # sorts an edge - replaced by function
    #S = lambda (x,y) : (x,y) if x<y else (y,x)
    
    # rs = Representative set of a vertex
    # for h in H, v in G is such that rs[h][v] == 1 if and only if v 
    # is a representant of h in G
    rs = p.new_variable(dim=2)
    
    for v in G:
        p.add_constraint(Sum([rs[h][v] for h in H]), max = 1)
    
    for h in H:
        for i in range(0,h):
            p.add_constraint(rs[h][i] == 0)

    # We ensure that the set of representatives of a
    # vertex h contains a tree, and thus is connected
    
    # edges represents the edges of the tree
    edges = p.new_variable(dim = 2)
    
    # there can be a edge for h between two vertices
    # only if those vertices represent h
    for u,v in G.edges(labels=None):
        for h in H:
            p.add_constraint(edges[h][S((u,v))] - rs[h][u], max = 0 )
            p.add_constraint(edges[h][S((u,v))] - rs[h][v], max = 0 )
    
    # The number of edges of the tree in h is exactly the cardinal
    # of its representative set minus 1
            
    for h in H:
        p.add_constraint(Sum([edges[h][S(e)] for e in G.edges(labels=None)])-Sum([rs[h][v] for v in G]), min=-1, max=-1)
    
    # a tree  has no cycle
    epsilon = 1/(5*Integer(G.order()))
    r_edges = p.new_variable(dim=2)
    
    for h in H:
        for u,v in G.edges(labels=None):
            p.add_constraint(r_edges[h][(u,v)] + r_edges[h][(v,u)] - edges[h][S((u,v))], min = 0)
                                                            
        for v in G:
            p.add_constraint(Sum([r_edges[h][(u,v)] for u in G.neighbors(v)]), max = 1-epsilon)
    
    # Once the representative sets are described, we must ensure
    # there are arcs corresponding to those of H between them
    h_edges = p.new_variable(dim=2)
    
    for h1, h2 in H.edges(labels=None):
    
        for v1, v2 in G.edges(labels=None):
    
            p.add_constraint(h_edges[(h1,h2)][S((v1,v2))] - rs[h2][v2], max = 0)
            p.add_constraint(h_edges[(h1,h2)][S((v1,v2))] - rs[h1][v1], max = 0)
            
            p.add_constraint(h_edges[(h2,h1)][S((v1,v2))] - rs[h1][v2], max = 0)
            p.add_constraint(h_edges[(h2,h1)][S((v1,v2))] - rs[h2][v1], max = 0)
        
        p.add_constraint(Sum([h_edges[(h1,h2)][S(e)] + h_edges[(h2,h1)][S(e)] for e in G.edges(labels=None) ]), min = 1)
    
    p.set_binary(rs)
    p.set_binary(edges)
    
    p.set_objective(None)
    
    try:
        p.solve(log=verbose)
    except MIPSolverException:
        raise ValueError("This graph has no minor isomorphic to H !")
    
    rs = p.get_values(rs)
    
    rs_dict = {}
    for h in H:
        rs_dict[h] = [v for v in G if rs[h][v]==1]
    
    return rs_dict