I am working on the Vertex Cover problem and trying to compare the performance of three different algorithms: Greedy (GRY), Pricing Algorithm (PA), and LP-based Rounding (LR). Specifically, I want to compare them over randomly generated instances with two parameters: n denoting the number of nodes taking values among {100, 120, 140, 160, 180, 200}, and p ∈ {0.1, 0.2, 0.3, 0.4, 0.5} representing the uniform probability of an edge existing among each pair of nodes.
For each given pair (n, p), I want to do the following sequentially:
- Generate M = 100 random graphs with each vertex cost uniformly sampled from [0, 1].
- For each randomly generated instance, implement the three algorithms and compute the ratio of the total cost on the vertex cover returned by the algorithm to that on the optimal, which can be obtained by solving an Integral Program (IP).
- Compute the average ratio over the M = 100 randomly generated instances for each of the three algorithms, denoted as τg(n, p), τp(n, p) and τl(n, p), which represent the performance of GRY, PA, and LR, respectively. However, I'm having trouble computing the ratios and comparing the performance of the algorithms.
How can I improve my implementation and comparison of the algorithms?
import random
import numpy as np
import matplotlib.pyplot as plt
import pulp
def generate_random_graph(n, p):
adj_matrix = [[0 for _ in range(n)] for _ in range(n)]
costs = [random.random() for _ in range(n)]
for i in range(n):
for j in range(i + 1, n):
if random.random() < p:
adj_matrix[i][j] = 1
adj_matrix[j][i] = 1
return adj_matrix, costs
def greedy_vc(adj_matrix, costs):
n = len(adj_matrix)
covered_edges = set()
vertex_cover = set()
while True:
uncovered_edges = [(i, j) for i in range(n) for j in range(i + 1, n) if adj_matrix[i][j] == 1 and (i, j) not in covered_edges]
if not uncovered_edges:
break
cost_ratio = [(i, costs[i] / sum(adj_matrix[i])) for i in range(n) if sum(adj_matrix[i]) > 0]
selected_vertex = min(cost_ratio, key=lambda x: x[1])[0]
vertex_cover.add(selected_vertex)
for j in range(n):
if adj_matrix[selected_vertex][j] == 1:
covered_edges.add((min(selected_vertex, j), max(selected_vertex, j)))
return vertex_cover
def pricing_algorithm_vc(adj_matrix, costs):
n = len(adj_matrix)
remaining_edges = {(i, j) for i in range(n) for j in range(i + 1, n) if adj_matrix[i][j] == 1}
vertex_cover = set()
while remaining_edges:
edge = remaining_edges.pop()
u, v = edge
if costs[u] < costs[v]:
vertex_cover.add(u)
for j in range(n):
if adj_matrix[u][j] == 1:
remaining_edges.discard((min(u, j), max(u, j)))
else:
vertex_cover.add(v)
for j in range(n):
if adj_matrix[v][j] == 1:
remaining_edges.discard((min(v, j), max(v, j)))
return vertex_cover
def lp_rounding_vc(adj_matrix, costs):
n = len(adj_matrix)
prob = pulp.LpProblem("VertexCoverLP", pulp.LpMinimize)
x_vars = [pulp.LpVariable(f"x_{i}", 0, 1, pulp.LpContinuous) for i in range(n)]
prob += pulp.lpSum(costs[i] * x_vars[i] for i in range(n))
for i in range(n):
for j in range(i + 1, n):
if adj_matrix[i][j] == 1:
prob += x_vars[i] + x_vars[j] >= 1
prob.solve()
vertex_cover = set(i for i in range(n) if x_vars[i].value() >= 0.5)
return vertex_cover
def optimal_vc(adj_matrix, costs):
n = len(adj_matrix)
prob = pulp.LpProblem("VertexCoverIP", pulp.LpMinimize)
x_vars = [pulp.LpVariable(f"x_{i}", 0, 1, pulp.LpInteger) for i in range(n)]
prob += pulp.lpSum(costs[i] * x_vars[i] for i in range(n))
for i in range(n):
for j in range(i + 1, n):
if adj_matrix[i][j] == 1:
prob += x_vars[i] + x_vars[j] >= 1
prob.solve()
vertex_cover = set(i for i in range(n) if x_vars[i].value() == 1)
return vertex_cover
def compare_algorithms(n_values, p_values):
champions = np.zeros((len(n_values), len(p_values)), dtype=int)
for n_idx, n in enumerate(n_values):
for p_idx, p in enumerate(p_values):
gry_ratios = []
pa_ratios = []
lr_ratios = []
for _ in range(10):
adj_matrix, costs = generate_random_graph(n, p)
gry_cover = greedy_vc(adj_matrix, costs)
pa_cover = pricing_algorithm_vc(adj_matrix, costs)
lr_cover = lp_rounding_vc(adj_matrix, costs)
opt_cover = optimal_vc(adj_matrix, costs)
if not opt_cover:
continue
gry_ratios.append(
sum(costs[v] for v in gry_cover) / sum(costs[v]
for v in opt_cover))
pa_ratios.append(
sum(costs[v] for v in pa_cover) / sum(costs[v]
for v in opt_cover))
lr_ratios.append(
sum(costs[v] for v in lr_cover) / sum(costs[v]
for v in opt_cover))
avg_gry = np.mean(gry_ratios)
avg_pa = np.mean(pa_ratios)
avg_lr = np.mean(lr_ratios)
min_ratio = min(avg_gry, avg_pa, avg_lr)
if min_ratio == avg_gry:
champions[n_idx][p_idx] = 1
elif min_ratio == avg_pa:
champions[n_idx][p_idx] = 2
else:
champions[n_idx][p_idx] = 3
return champions
def plot_champions(champions, n_values, p_values):
colors = {1: 'red', 2: 'blue', 3: 'black'}
for i, n in enumerate(n_values):
for j, p in enumerate(p_values):
plt.scatter(n, p, color=colors[champions[i][j]], marker='s')
plt.xlabel('n')
plt.ylabel('p')
plt.show()
n_values = [2,3]
p_values = [0.1, 0.2]
champions = compare_algorithms(n_values, p_values)
plot_champions(champions, n_values, p_values)
