# Solving vertex cover with a quantum annealer

The problem of vertex cover is, given an undirected graph $G = (V, E)$, colour the smallest amount of vertices such that each edge $e \in E$ is connected to a coloured vertex.

This notebooks works through the process of creating a random graph, translating to an optimization problem, and eventually finding the ground state using a quantum annealer.

### Graph setup

The first thing we will do is create an instance of the problem, by constructing a small, random undirected graph. We are going to use the `networkx` package, which should already be installed if you have installed if you are using Anaconda.


```python
import dimod
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
```


```python
n_vertices = 5
n_edges = 6

small_graph = nx.gnm_random_graph(n_vertices, n_edges)
```


```python
nx.draw(small_graph, with_labels=True)
```

### Constructing the Hamiltonian

I showed in class that the objective function for vertex cover looks like this:
\begin{equation}
 \sum_{(u,v) \in E} (1 - x_u) (1 - x_v) + \gamma \sum_{v \in V} x_v
\end{equation}
We want to find an assignment of the $x_u$ of 1 (coloured) or 0 (uncoloured) that _minimizes_ this function. The first sum tries to force us to choose an assignment that makes sure every edge gets attached to a coloured vertex. The second sum is essentially just counting the number of coloured vertices.

**Task**: Expand out the QUBO above to see how you can convert it to a more 'traditional' looking QUBO:
\begin{equation}
 \sum_{(u,v) \in E} x_u x_v + \sum_{v \in V} (\gamma - \text{deg}(x_v)) x_v
\end{equation}
where deg($x_v$) indicates the degree of vertex $x_v$ in the graph.


```python
γ = 0.8
Q = {x : 1 for x in small_graph.edges()}
r = {x : (γ - small_graph.degree[x]) for x in small_graph.nodes}
```

Let's convert it to the appropriate data structure, and solve using the exact solver. 


```python
bqm = dimod.BinaryQuadraticModel(r, Q, 0, dimod.BINARY)
response = dimod.ExactSolver().sample(bqm)
print(f"Sample energy = {next(response.data(['energy']))[0]}")
```

Let's print the graph with proper colours included


```python
colour_assignments = next(response.data(['sample']))[0]
colours = ['grey' if colour_assignments[x] == 0 else 'red' for x in range(len(colour_assignments))]

nx.draw(small_graph, with_labels=True, node_color=colours)
```

### Scaling up...

That one was easy enough to solve by hand. Let's try a much larger instance...


```python
n_vertices = 20
n_edges = 60

large_graph = nx.gnm_random_graph(n_vertices, n_edges)
nx.draw(large_graph, with_labels=True)
```


```python
# Create h, J and put it into the exact solver
γ = 0.8
Q = {x : 1 for x in large_graph.edges()}
r = {x : (γ - large_graph.degree[x]) for x in large_graph.nodes}

bqm = dimod.BinaryQuadraticModel(r, Q, 0, dimod.BINARY)
response = dimod.ExactSolver().sample(bqm)
print(f"Sample energy = {next(response.data(['energy']))[0]}")
      
colour_assignments = next(response.data(['sample']))[0]
colours = ['grey' if colour_assignments[x] == 0 else 'red' for x in range(len(colour_assignments))]

nx.draw(large_graph, with_labels=True, node_color=colours)
print(f"Coloured {list(colour_assignments.values()).count(1)}/{n_vertices} vertices.")
```

### Running on the D-Wave

You'll only be able to run the next few cells if you have D-Wave access. We will send the same graph as before to the D-Wave QPU and see what kind of results we get back!


```python
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite

sampler = EmbeddingComposite(DWaveSampler())
```


```python
ising_conversion = bqm.to_ising()
h, J = ising_conversion[0], ising_conversion[1]
response = sampler.sample_ising(h, J, num_reads = 1000)
```


```python
best_solution =np.sort(response.record, order='energy')[0]
```


```python
print(f"Sample energy = {best_solution['energy']}")
      
colour_assignments_qpu = {x : best_solution['sample'][x] for x in range(n_vertices)}
for x in range(n_vertices):
      if colour_assignments_qpu[x] == -1:
          colour_assignments_qpu[x] = 0
colours = ['grey' if colour_assignments_qpu[x] == 0 else 'red' for x in range(len(colour_assignments_qpu))]

nx.draw(large_graph, with_labels=True, node_color=colours)
print(f"Coloured {list(colour_assignments_qpu.values()).count(1)}/{n_vertices} vertices.")
```


```python
print("Node\tExact\tQPU")
for x in range(n_vertices):
    print(f"{x}\t{colour_assignments[x]}\t{colour_assignments_qpu[x]}")
```

Here is a scatter plot of all the different energies we got out, against the number of times each solution occurred. 


```python
plt.scatter(response.record['energy'], response.record['num_occurrences'])
```


```python
response.record['num_occurrences']
```

269e15dbd63116b4c521a64b0e0ac8fff41fca84
ipynb
Jupyter Notebook
04-annealing-applications/Vertex-Cover.ipynb
glassnotes/Intro-QC-TRIUMF
1b54a7fd9795ff5a8541c9a99fdb75365cbd113a
MIT
2019-05-09T17:40:20.000Z
2021-12-15T12:23:17.000Z
04-annealing-applications/Vertex-Cover.ipynb
glassnotes/Intro-QC-TRIUMF
1b54a7fd9795ff5a8541c9a99fdb75365cbd113a
MIT
2021-09-29T07:34:09.000Z
2021-09-29T21:01:29.000Z
04-annealing-applications/Vertex-Cover.ipynb
glassnotes/Intro-QC-TRIUMF
1b54a7fd9795ff5a8541c9a99fdb75365cbd113a
MIT
2019-05-09T18:45:49.000Z
2021-12-15T12:23:21.000Z
Qwen/Qwen-72B
1. YES
2. YES