Merge pull request #1002 from qiboteam/mvc2.0

fixing previous mvc branch

doc/source/code-examples/applications.rst

@@ -93,6 +93,7 @@ Combinatorics
     :maxdepth: 1
 
     tutorials/qap/README.md
+    tutorials/mvc/README.md
 
 
 Applications by algorithm

doc/source/code-examples/tutorials/mvc/README.md

@@ -0,0 +1 @@
+../../../../../examples/mvc/README.md

examples/README.md

@@ -27,6 +27,7 @@ physics problems.
 - [Quantum k-medians clustering](qclustering/README.md)
 - [Determining probability density functions with adiabatic quantum computing](adiabatic_qml/adiabatic-qml.ipynb)
 - [Quadratic assignment problem (QAP)](qap/README.md)
+- [Minimum Vertex Cover (MVC)](mvc/README.md)
 
 In the `benchmarks` folder we have included examples concerning:
 - A generic benchmark script for multiple circuits (`benchmarks/main.py`)

examples/mvc/README.md

@@ -0,0 +1,228 @@
+# Minimum Vertex Cover (MVC)
+
+Code at: [https://github.com/qiboteam/qibo/tree/master/examples/mvc](https://github.com/qiboteam/qibo/tree/master/examples/mvc)
+
+Given an undirected graph with a set of nodes V and edges E, A vertex cover of a graph is a set of nodes that includes at
+least one endpoint of every edge of the graph. The Minimum Vertex Cover (MVC) problem is an optimisation problem
+that find smallest vertex cover of a given graph.
+
+## Load graph from csv
+
+
+```python
+import networkx as nx
+import csv
+import matplotlib.pyplot as plt
+
+from mvc import qubo_mvc, qubo_mvc_penalty, mvc_feasibility, mvc_easy_fix, mvc_energy
+
+
+def load_csv(filename:str):
+    """ Load graph from csv file
+    """
+
+    with open(filename, 'r', newline='') as f:
+        reader = csv.reader(f)
+        data = [[int(row[0]), int(row[1]), float(row[2])] for row in reader]
+
+    nodes = [k0 for k0, k1, v in data if k0 == k1]
+    edges = [[k0, k1] for k0, k1, v in data if k0 != k1]
+
+    g = nx.Graph()
+    g.add_nodes_from(nodes)
+    g.add_edges_from(edges)
+
+    node_weights = {k0: {'weight': v} for k0, k1, v in data if k0 == k1}
+    edge_weights = {(k0, k1): {'weight': v} for k0, k1, v in data if k0 != k1}
+
+    nx.set_edge_attributes(g, edge_weights)
+    nx.set_node_attributes(g, node_weights)
+
+    return g
+```
+
+csv Format:
+
+0,0,0.4696822179283675
+1,1,0.6917392308135916
+2,2,0.7130420958314817
+3,3,0.49342677332980056
+4,4,0.49661600571433673
+5,5,0.55601135361347
+6,6,0.2831095711244086
+7,7,0.18737823232636885
+0,1,0.98602725407379
+1,2,0.935853268681996
+2,3,0.23080434023289664
+4,5,0.25521731195623476
+5,6,0.37096743935296606
+6,7,0.04408120693490547
+0,2,0.5619060764177991
+4,7,0.48402095290884917
+1,6,0.7106584134580318
+
+
+
+```python
+g = load_csv('mvc.csv')
+```
+
+
+```python
+pos = nx.spring_layout(g, seed=1234)
+nx.draw(g, pos)
+```
+
+
+
+![png](README_files/README_7_0.png)
+
+
+
+
+```python
+print(f'Number of node is {len(g.nodes)}')
+print(f'Number of edge is {len(g.edges)}')
+```
+
+    Number of node is 8
+    Number of edge is 9
+
+
+## QUBO formulation
+
+### Estimate penalty
+
+
+```python
+penalty = qubo_mvc_penalty(g)
+print(f'The penalty is {penalty}')
+```
+
+    The penalty is 1.6548482220717595
+
+
+### Formulate QUBO (weighted minimum vertex cover)
+
+
+```python
+linear, quadratic = qubo_mvc(g, penalty=penalty)
+print(f'Number of linear terms: {len(linear)}')
+print(f'Number of quadratic terms: {len(quadratic)}\n')
+```
+
+    Number of linear terms: 8
+    Number of quadratic terms: 9
+
+
+
+### Generate random solutions
+
+
+```python
+import numpy as np
+rng = np.random.default_rng(seed=1234)
+random_solution = {i: rng.integers(2) for i in g.nodes}
+print(f'The random solution is {random_solution}\n')
+```
+
+    The random solution is {0: 1, 1: 1, 2: 1, 3: 0, 4: 0, 5: 1, 6: 0, 7: 0}
+
+
+
+### Check feasibility
+
+
+```python
+feasibility = mvc_feasibility(g, random_solution)
+print(f'The feasibility is {feasibility}\n')
+```
+
+    The feasibility is False
+
+
+
+### Fix broken constraints
+
+
+```python
+fixed_solution = mvc_easy_fix(g, random_solution)
+print(f'After fixation, the solution is {fixed_solution}\n')
+```
+
+    After fixation, the solution is {0: 1, 1: 1, 2: 1, 3: 0, 4: 1, 5: 1, 6: 1, 7: 0}
+
+
+
+
+```python
+feasibility = mvc_feasibility(g, fixed_solution)
+print(f'The feasibility is {feasibility}\n')
+```
+
+    The feasibility is True
+
+
+
+### Visualisation
+
+
+```python
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12.8,4.8))
+colour_b4fix = ['C1' if random_solution[n] else 'C0' for n in g.nodes]
+colour_fixed = ['C1' if fixed_solution[n] else 'C0' for n in g.nodes]
+nx.draw(g, pos, node_color=colour_b4fix, ax=ax0)
+ax0.set_title(f'Vertex cover in orange before fix')
+nx.draw(g, pos, node_color=colour_fixed, ax=ax1)
+ax1.set_title(f'Vertex cover in orange after fix')
+```
+
+
+
+
+    Text(0.5, 1.0, 'Vertex cover in orange after fix')
+
+
+
+
+
+![png](README_files/README_22_1.png)
+
+
+
+### Calculate energy
+
+
+```python
+energy = mvc_energy(g, fixed_solution)
+print(f'The energy is {energy}')
+```
+
+    The energy is 3.2102004750256556
+
+
+##  Hamiltonian of the MVC problem
+
+
+```python
+ham = hamiltonian_mvc(g, penalty=penalty, dense=True)
+```
+
+    [Qibo 0.1.8|INFO|2022-11-01 10:26:19]: Using numpy backend on /CPU:0
+
+
+## Solve the hamiltonian with QAOA
+
+
+```python
+from qibo import models, hamiltonians
+
+# Create QAOA model
+qaoa = models.QAOA(ham)
+
+# Optimize starting from a random guess for the variational parameters
+initial_parameters = 0.01 * np.random.uniform(0,1,2)
+best_energy, final_parameters, extra = qaoa.minimize(initial_parameters, method="BFGS")
+```
+
+    [Qibo 0.1.8|WARNING|2022-10-31 14:14:37]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.

examples/mvc/main.py

@@ -0,0 +1,81 @@
+"""Minimum Vertex Cover"""
+import argparse
+import csv
+
+import networkx as nx
+from mvc import (
+    hamiltonian_mvc,
+    mvc_easy_fix,
+    mvc_energy,
+    mvc_feasibility,
+    qubo_mvc,
+    qubo_mvc_penalty,
+)
+from qubo_utils import binary2spin, spin2QiboHamiltonian
+
+from qibo import callbacks, hamiltonians, models
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--filename", default="./mvc.csv", type=str)
+
+
+def load_csv(filename: str):
+    """Load graph from csv file"""
+    with open(filename, newline="") as f:
+        reader = csv.reader(f)
+        data = [[int(row[0]), int(row[1]), float(row[2])] for row in reader]
+
+    nodes = [k0 for k0, k1, v in data if k0 == k1]
+    edges = [[k0, k1] for k0, k1, v in data if k0 != k1]
+
+    g = nx.Graph()
+    g.add_nodes_from(nodes)
+    g.add_edges_from(edges)
+
+    node_weights = {k0: {"weight": v} for k0, k1, v in data if k0 == k1}
+    edge_weights = {(k0, k1): {"weight": v} for k0, k1, v in data if k0 != k1}
+
+    nx.set_edge_attributes(g, edge_weights)
+    nx.set_node_attributes(g, node_weights)
+
+    return g
+
+
+def main(filename: str = "./mvc.csv"):
+    print(f"Load a graph from {filename} and make it a QUBO")
+    g = load_csv(filename)
+    penalty = qubo_mvc_penalty(g)
+    linear, quadratic = qubo_mvc(g, penalty=penalty)
+
+    print("A random solution with seed 1234 must be infeasible")
+    import numpy as np
+
+    rng = np.random.default_rng(seed=1234)
+    random_solution = {i: rng.integers(2) for i in g.nodes}
+    feasibility = mvc_feasibility(g, random_solution)
+    assert not feasibility, "The random solution should be infeasible."
+
+    print("Make a naive fix to the violation of the constraint")
+    fixed_solution = mvc_easy_fix(g, random_solution)
+    feasibility = mvc_feasibility(g, fixed_solution)
+    assert feasibility, "The fixed solution should be feasible."
+
+    print("Calculate the energy of the solution")
+    energy = mvc_energy(g, fixed_solution)
+
+    print("Construct a hamiltonian based on the QUBO representation")
+    h, J, _ = binary2spin(linear, quadratic)
+    h = {k: -v for k, v in h.items()}
+    J = {k: -v for k, v in J.items()}
+    ham = spin2QiboHamiltonian(h, J)
+
+    print("Construct a hamiltonian directly from a networkx graph")
+    ham = hamiltonian_mvc(g)
+
+    print("done.")
+
+
+if __name__ == "__main__":
+    # by defualt, test on the mvc.csv in the same directory
+    args = parser.parse_args()
+    main(args.filename)

examples/mvc/mvc.csv

@@ -0,0 +1,17 @@
+0,0,0.4696822179283675
+1,1,0.6917392308135916
+2,2,0.7130420958314817
+3,3,0.49342677332980056
+4,4,0.49661600571433673
+5,5,0.55601135361347
+6,6,0.2831095711244086
+7,7,0.18737823232636885
+0,1,0.98602725407379
+1,2,0.935853268681996
+2,3,0.23080434023289664
+4,5,0.25521731195623476
+5,6,0.37096743935296606
+6,7,0.04408120693490547
+0,2,0.5619060764177991
+4,7,0.48402095290884917
+1,6,0.7106584134580318