Maximum weighted clique

maxclique.py

"""
Maximum weighted clique for a graph
"""


import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import networkx as nx
import cplex
from cplex.callbacks import LazyConstraintCallback


def greedy_expand(G, init_set):

    R = G.copy()
    neigh = [n for i in init_set for n in R.neighbors(i)]
    R.remove_nodes_from(init_set)
    R.remove_nodes_from(neigh)
    if R.number_of_nodes() != 0:
        x = get_min_degree_vertex(R)
    while R.number_of_nodes() != 0:
        neigh2 = [m for m in R.neighbors(x)]
        R.remove_node(x)
        init_set.append(x)
        R.remove_nodes_from(neigh2)
        if R.number_of_nodes() != 0:
            x = get_min_degree_vertex(R)

    return init_set


def greedy_init(G):
    """
    https://kmutya.github.io/maxclique/
    """
    n = G.number_of_nodes()  # Storing total number of nodes in 'n'
    max_ind_sets = []  # initializing a list that will store maximum independent sets
    for j in G.nodes:
        R = G.copy()  # Storing a copy of the graph as a residual
        neigh = [n for n in R.neighbors(j)]  # Catch all the neighbours of j
        R.remove_node(j)  # removing the node we start from
        iset = [j]
        R.remove_nodes_from(neigh)  # Removing the neighbours of j
        if R.number_of_nodes() != 0:
            x = get_min_degree_vertex(R)
        while R.number_of_nodes() != 0:
            neigh2 = [m for m in R.neighbors(x)]
            R.remove_node(x)
            iset.append(x)
            R.remove_nodes_from(neigh2)
            if R.number_of_nodes() != 0:
                x = get_min_degree_vertex(R)

        max_ind_sets.append(iset)

    return(max_ind_sets)


def get_min_degree_vertex(Residual_graph):
    '''Takes in the residual graph R and returns the node with the lowest degree'''
    degrees = [val for (node, val) in Residual_graph.degree()]
    node = [node for (node, val) in Residual_graph.degree()]
    node_degree = dict(zip(node, degrees))
    return (min(node_degree, key=node_degree.get))


G = nx.generators.ring_of_cliques(6, 3)

nx.draw(G, with_labels=True)
plt.savefig("tmp.png")
mis = greedy_init(G)

prob = cplex.Cplex()

numvar = len(G.nodes)


def func(x): return "x"+str(x)


class LazyCallback(LazyConstraintCallback):
    """Lazy constraint callback to generate maximum independent sets on the fly.

    There are too many such constraints to make them all available to 
    CPLEX right away - and in any case, very few of them are valid.

    So generate them on the fly.
    """

    # Callback constructor.
    #
    # Any needed fields are set externally after registering the callback.
    def __init__(self, env):
        super().__init__(env)

    def __call__(self):

        values = self.get_values(self.names)

        # Any node with non-zero value is considered as part of the set
        curr_solution = [int(name[1:]) for name, val in zip(self.names, values) if val >= 0.001]
        print("Current solution: ", curr_solution)

        # Look to generate all independent sets that would cut off the (fractional)
        # value. To do this, first induce a subgraph - and for each node, built a
        #
        subG = self.G.subgraph(curr_solution)
        sub_ind_set = greedy_init(subG)
        max_ind_sets = [greedy_expand(self.G, sset) for sset in sub_ind_set]

        for iset in max_ind_sets:

            con_vars = [func(i) for i in iset]
            coeffs = [1.0] * len(con_vars)
            lhs = cplex.SparsePair(con_vars, coeffs)
            self.add(constraint=lhs, rhs=1.0, sense='L')


names = list(map(func, G.nodes))
var_type = [prob.variables.type.continuous] * numvar
prob.variables.add(names=names,
                   lb=[0.0] * numvar,
                   ub=[1.0] * numvar,
                   types=var_type)
prob.objective.set_sense(prob.objective.sense.maximize)
prob.objective.set_linear([(n, 1.0) for n in names])
lhs = []
for iset in mis:
    con_vars = [func(i) for i in iset]
    coeffs = [1.0] * len(con_vars)
    lhs.append(cplex.SparsePair(con_vars, coeffs))
prob.linear_constraints.add(lin_expr=lhs,
                            rhs=[1.0] * len(lhs),
                            senses=['L'] * len(lhs))
print("Constraint: Maximum independent set. ({} constraints)".format(len(mis)))

# register callbacks to generate additional independent sets on the fly
lazycb = prob.register_callback(LazyCallback)
lazycb.names = names
lazycb.G = G

# prob.write("test.lp")
prob.solve()
status = prob.solution.status[prob.solution.get_status()]
print("Status:{}".format(status))

if prob.solution.get_status() in [101, 105, 107, 111, 113]:
    # Optimal/feasible

    # Get the solution
    print("Solution value: ")
    print(prob.solution.get_objective_value())

    # get the configuration
    x_res = prob.solution.get_values(names)
    for x_name, val in zip(names, x_res):
        if val > 0:
            print(x_name, val)

Constraint 2: All maximum independent set. (18 constraints)
Version identifier: 12.10.0.0 | 2019-11-26 | 843d4de
CPXPARAM_Read_DataCheck 1
Warning: Control callbacks may disable some MIP features.
Lazy constraint(s) or lazy constraint/branch callback is present.
Disabling dual reductions (CPX_PARAM_REDUCE) in presolve.
Disabling non-linear reductions (CPX_PARAM_PRELINEAR) in presolve.
Disabling presolve reductions that prevent crushing forms.
Disabling repeat represolve because of lazy constraint/incumbent callback.
Tried aggregator 1 time.
MIP Presolve eliminated 6 rows and 0 columns.
Reduced MIP has 12 rows, 18 columns, and 72 nonzeros.
Reduced MIP has 0 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.02 ticks)
MIP emphasis: balance optimality and feasibility.
MIP search method: traditional branch-and-cut.
Parallel mode: none, using 1 thread.
Root relaxation solution time = 0.00 sec. (0.01 ticks)
Current solution: [5, 7, 8, 14, 16, 17]
Current solution: [4, 10, 11, 13]
Current solution: [0, 6, 8, 10, 11]
Current solution: [4, 6, 8, 10, 11, 16]
Current solution: [0, 2, 4, 5, 6, 8]
Current solution: [6, 7, 8]
Current solution: [6, 7, 8]

    Nodes                                         Cuts/

Node Left Objective IInf Best Integer Best Bound ItCnt Gap Variable B NodeID Parent Depth

• 0     0      integral     0        3.0000        6.0000        8  100.00%                        0             0
  

Elapsed time = 0.02 sec. (0.16 ticks, tree = 0.00 MB, solutions = 1)

User cuts applied: 17

Root node processing (before b&c):
Real time = 0.02 sec. (0.16 ticks)
Sequential b&c:
Real time = 0.00 sec. (0.00 ticks)
------------
Total (root+branch&cut) = 0.02 sec. (0.16 ticks)
Status:MIP_optimal
Solution value:
3.0
x6 1.0
x7 1.0
x8 1.0