annealing_test_mds.py

# Description：simulated annealing for 2-hop minimum dominating set
import networkx as nx
import random
import math
from collections import deque
import matplotlib.pyplot as plt
import time

def create_random_digraph(n_nodes, edge_probability=0.2):
    G = nx.DiGraph()
    G.add_nodes_from(range(n_nodes))
    
    for i in range(n_nodes):
        for j in range(n_nodes):
            if i != j and random.random() < edge_probability:
                G.add_edge(i, j)
    
    while not nx.is_strongly_connected(G):
        i, j = random.sample(range(n_nodes), 2)
        if not G.has_edge(i, j):
            G.add_edge(i, j)
    
    return G

class TwoHopMDS:
    def __init__(self, G):
        self.G = G
        self.n = G.number_of_nodes()
        self.two_hop_reachable = self._compute_two_hop_reachable()
    
    def _compute_two_hop_reachable(self):
        reachable = {}
        for node in self.G.nodes():
            reached = {node}
            
            first_hop = set(self.G.successors(node))
            reached.update(first_hop)
            
            for successor in first_hop:
                reached.update(self.G.successors(successor))
            
            reachable[node] = reached
        return reachable
    
    def is_two_hop_dominating(self, nodes):
        if not nodes:
            return False
        
        covered = set()
        for node in nodes:
            covered.update(self.two_hop_reachable[node])
        
        return len(covered) == self.n
    
    def get_initial_solution(self):
        solution = set()
        uncovered = set(self.G.nodes())
        
        while uncovered:
            max_coverage = -1
            best_node = None
            
            for node in self.G.nodes():
                if node not in solution:
                    coverage = len(self.two_hop_reachable[node] & uncovered)
                    if coverage > max_coverage:
                        max_coverage = coverage
                        best_node = node
            
            if best_node is None:
                break
                
            solution.add(best_node)
            uncovered -= self.two_hop_reachable[best_node]
        
        return solution
    
    def simulated_annealing(self, initial_temp=100, cooling_rate=0.95, iterations_per_temp=100):
        def get_neighbor(current):
            neighbor = current.copy()
            operation = random.random()
            
            if operation < 0.4 and len(neighbor) > 1:
                node_to_remove = random.choice(list(neighbor))
                neighbor.remove(node_to_remove)
            elif operation < 0.8:
                available_nodes = set(self.G.nodes()) - neighbor
                if available_nodes:
                    node_to_add = random.choice(list(available_nodes))
                    neighbor.add(node_to_add)
            else:
                if neighbor and len(self.G.nodes()) > len(neighbor):
                    node_to_remove = random.choice(list(neighbor))
                    available_nodes = set(self.G.nodes()) - neighbor
                    node_to_add = random.choice(list(available_nodes))
                    neighbor.remove(node_to_remove)
                    neighbor.add(node_to_add)
            
            return neighbor

        current_solution = self.get_initial_solution()
        best_solution = current_solution.copy()
        current_temp = initial_temp
        no_improvement = 0
        
        while current_temp > 0.1 and no_improvement < 100:
            improved = False
            
            for _ in range(iterations_per_temp):
                new_solution = get_neighbor(current_solution)
                
                if self.is_two_hop_dominating(new_solution):
                    delta = len(new_solution) - len(current_solution)
                    
                    if delta < 0:
                        current_solution = new_solution
                        if len(new_solution) < len(best_solution):
                            best_solution = new_solution.copy()
                            improved = True
                    else:
                        acceptance_probability = math.exp(-delta / current_temp)
                        if random.random() < acceptance_probability:
                            current_solution = new_solution
            
            current_temp *= cooling_rate
            if not improved:
                no_improvement += 1
            else:
                no_improvement = 0
        
        return best_solution

def visualize_reachability(G, node, pos=None):
    if pos is None:
        pos = nx.spring_layout(G)
    
    plt.figure(figsize=(12, 8))
    
    first_hop = set(G.successors(node))
    second_hop = set()
    for n in first_hop:
        second_hop.update(G.successors(n))
    second_hop -= first_hop
    second_hop.discard(node)
    
    nx.draw_networkx_edges(G, pos, edge_color='gray', alpha=0.2)
    
    edges_first_hop = [(node, n) for n in first_hop]
    nx.draw_networkx_edges(G, pos, edgelist=edges_first_hop, 
                          edge_color='blue', width=2)
    
    edges_second_hop = []
    for n1 in first_hop:
        for n2 in G.successors(n1):
            if n2 in second_hop:
                edges_second_hop.append((n1, n2))
    nx.draw_networkx_edges(G, pos, edgelist=edges_second_hop, 
                          edge_color='green', width=2)
    
    other_nodes = set(G.nodes()) - {node} - first_hop - second_hop
    nx.draw_networkx_nodes(G, pos, nodelist=[node], node_color='red', 
                          node_size=500, label='Source')
    nx.draw_networkx_nodes(G, pos, nodelist=list(first_hop), 
                          node_color='blue', node_size=500, label='First hop')
    nx.draw_networkx_nodes(G, pos, nodelist=list(second_hop), 
                          node_color='green', node_size=500, label='Second hop')
    nx.draw_networkx_nodes(G, pos, nodelist=list(other_nodes), 
                          node_color='gray', node_size=500, label='Unreachable')
    
    nx.draw_networkx_labels(G, pos)
    
    plt.title(f"2-Hop Reachability from Node {node}")
    plt.legend()
    plt.axis('off')
    plt.show()

def test_with_random_digraph(n_nodes=15, edge_prob=0.2):
    G = create_random_digraph(n_nodes, edge_prob)
    
    start_time = time.time()
    
    solver = TwoHopMDS(G)
    solution = solver.simulated_annealing(
        initial_temp=100,
        cooling_rate=0.95,
        iterations_per_temp=50
    )
    
    runtime = time.time() - start_time
    
    print(f"graph info:")
    print(f"node number: {G.number_of_nodes()}")
    print(f"edge number: {G.number_of_edges()}")
    print(f"\nsolve result:")
    print(f"2-hop dominating set: {sorted(list(solution))}")
    print(f"set size: {len(solution)}")
    print(f"is valid solution: {solver.is_two_hop_dominating(solution)}")
    print(f"run time: {runtime:.2f} seconds")
    
    # visualize the graph and solution
    pos = nx.spring_layout(G, k=1, iterations=50)
    
    plt.figure(figsize=(12, 8))
    nx.draw_networkx_edges(G, pos, edge_color='gray', 
                          arrowsize=20, width=1, alpha=0.5)
    
    non_dominating = set(G.nodes()) - solution
    nx.draw_networkx_nodes(G, pos, nodelist=list(non_dominating),
                          node_color='lightblue', node_size=500, alpha=0.6)
    
    nx.draw_networkx_nodes(G, pos, nodelist=list(solution),
                          node_color='red', node_size=500)
    
    nx.draw_networkx_labels(G, pos)
    
    plt.title("2-Hop Dominating Set in Random Directed Graph")
    plt.axis('off')
    plt.show()
    
    for node in solution:
        visualize_reachability(G, node, pos)
    
    return G, solution

# only for testing
if __name__ == "__main__":
    random.seed(42)  
    G, solution = test_with_random_digraph(500,0.01)