diff --git a/doc/source/code-examples/applications.rst b/doc/source/code-examples/applications.rst index 254b164e6e..ede7bb33bc 100644 --- a/doc/source/code-examples/applications.rst +++ b/doc/source/code-examples/applications.rst @@ -86,6 +86,15 @@ Quantum Machine Learning tutorials/qclustering/README.md tutorials/adiabatic_qml/adiabatic-qml.ipynb +Combinatorics +^^^^^^^^^^^^^ + +.. toctree:: + :maxdepth: 1 + + tutorials/qap/README.md + + Applications by algorithm ------------------------- diff --git a/doc/source/code-examples/tutorials/qap/README.md b/doc/source/code-examples/tutorials/qap/README.md new file mode 120000 index 0000000000..047aaaa93c --- /dev/null +++ b/doc/source/code-examples/tutorials/qap/README.md @@ -0,0 +1 @@ +../../../../../examples/qap/README.md \ No newline at end of file diff --git a/examples/README.md b/examples/README.md index 4c32c815a6..17ca621ed2 100644 --- a/examples/README.md +++ b/examples/README.md @@ -26,6 +26,7 @@ physics problems. - [Quantum anomaly detection](anomaly_detection/README.md) - [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) In the `benchmarks` folder we have included examples concerning: - A generic benchmark script for multiple circuits (`benchmarks/main.py`) diff --git a/examples/qap/README.md b/examples/qap/README.md new file mode 100644 index 0000000000..d29e621de0 --- /dev/null +++ b/examples/qap/README.md @@ -0,0 +1,177 @@ +# Quadratic assignment problem (QAP) + +Code at: [https://github.com/qiboteam/qibo/tree/master/examples/qap](https://github.com/qiboteam/qibo/tree/master/examples/qap) + +The quadratic assignment problem (QAP) is an important combinatorial optimization problems that was first introduced by Koopmans and Beckmann. The objective of the problem is to assign a set of facilities to a set of locations in such a way as to minimize the total assignment cost. The assignment cost for a pair of facilities is a function of the flow between the facilities and the distance between the locations of the facilities. + +```python +import numpy as np + +from qap import qubo_qap, qubo_qap_penalty, qubo_qap_feasibility, qubo_qap_energy, hamiltonian_qap + +def load_qap(filename): + """Load qap problem from a file + + The file format is compatible with the one used in QAPLIB + + """ + + with open(filename, 'r') as fh: + n = int(fh.readline()) + + numbers = [float(n) for n in fh.read().split()] + + data = np.asarray(numbers).reshape(2, n, n) + f = data[1] + d = data[0] + + i = range(len(f)) + f[i, i] = 0 + d[i, i] = 0 + + return f, d +``` + +## Load QAP problem from a file + + +```python +F, D = load_qap('tiny04a.dat') +print(f'The QAP instance is:') +print(F) +print(D) +``` + + The QAP instance is: + [[0. 0.29541331 0.68442855 0.19882279] + [0.29541331 0. 0.61649225 0.16210679] + [0.68442855 0.61649225 0. 0.73052088] + [0.19882279 0.16210679 0.73052088 0. ]] + [[0. 0.77969778 0.43045022 0.43294055] + [0.77969778 0. 0.1920096 0.58829618] + [0.43045022 0.1920096 0. 0.47901122] + [0.43294055 0.58829618 0.47901122 0. ]] + + +## Calculate the penalty + + +```python +penalty = qubo_qap_penalty((F, D)) +print(f'The penalty is {penalty}') +``` + + The penalty is 2.2783420340595995 + + +## Formulate the QUBO + + +```python +linear, quadratic, offset = qubo_qap((F, D), penalty=penalty) +print(f'linear: {linear}') +print() +print(f'quadratic: {quadratic}') +print() +print(f'offset: {offset}\n') +``` + + linear: {0: -4.556684, 1: -4.556684, 2: -4.556684, 3: -4.556684, 4: -4.556684, 5: -4.556684, 6: -4.556684, 7: -4.556684, 8: -4.556684, 9: -4.556684, 10: -4.556684, 11: -4.556684, 12: -4.556684, 13: -4.556684, 14: -4.556684, 15: -4.556684} + + quadratic: {(1, 0): 2.278342, (2, 0): 2.278342, (3, 0): 2.278342, (4, 0): 2.278342, (5, 0): 0.2303331, (6, 0): 0.12716073, (7, 0): 0.1278964, (8, 0): 2.278342, (9, 0): 0.5336474, (10, 0): 0.2946124, (11, 0): 0.29631686, (12, 0): 2.278342, (13, 0): 0.15502168, (14, 0): 0.085583314, (15, 0): 0.08607845, (2, 1): 2.278342, (3, 1): 2.278342, (4, 1): 0.2303331, (5, 1): 2.278342, (6, 1): 0.05672219, (7, 1): 0.17379051, (8, 1): 0.5336474, (9, 1): 2.278342, (10, 1): 0.13141686, (11, 1): 0.4026467, (12, 1): 0.15502168, (13, 1): 2.278342, (14, 1): 0.038175885, (15, 1): 0.11696669, (3, 2): 2.278342, (4, 2): 0.12716073, (5, 2): 0.05672219, (6, 2): 2.278342, (7, 2): 0.14150628, (8, 2): 0.2946124, (9, 2): 0.13141686, (10, 2): 2.278342, (11, 2): 0.32784894, (12, 2): 0.085583314, (13, 2): 0.038175885, (14, 2): 2.278342, (15, 2): 0.09523835, (4, 3): 0.1278964, (5, 3): 0.17379051, (6, 3): 0.14150628, (7, 3): 2.278342, (8, 3): 0.29631686, (9, 3): 0.4026467, (10, 3): 0.32784894, (11, 3): 2.278342, (12, 3): 0.08607845, (13, 3): 0.11696669, (14, 3): 0.09523835, (15, 3): 2.278342, (5, 4): 2.278342, (6, 4): 2.278342, (7, 4): 2.278342, (8, 4): 2.278342, (9, 4): 0.48067763, (10, 4): 0.2653692, (11, 4): 0.2669045, (12, 4): 2.278342, (13, 4): 0.1263943, (14, 4): 0.069778904, (15, 4): 0.0701826, (6, 5): 2.278342, (7, 5): 2.278342, (8, 5): 0.48067763, (9, 5): 2.278342, (10, 5): 0.11837243, (11, 5): 0.36268005, (12, 5): 0.1263943, (13, 5): 2.278342, (14, 5): 0.03112606, (15, 5): 0.095366806, (7, 6): 2.278342, (8, 6): 0.2653692, (9, 6): 0.11837243, (10, 6): 2.278342, (11, 6): 0.2953067, (12, 6): 0.069778904, (13, 6): 0.03112606, (14, 6): 2.278342, (15, 6): 0.07765097, (8, 7): 0.2669045, (9, 7): 0.36268005, (10, 7): 0.2953067, (11, 7): 2.278342, (12, 7): 0.0701826, (13, 7): 0.095366806, (14, 7): 0.07765097, (15, 7): 2.278342, (9, 8): 2.278342, (10, 8): 2.278342, (11, 8): 2.278342, (12, 8): 2.278342, (13, 8): 0.5695855, (14, 8): 0.31445286, (15, 8): 0.3162721, (10, 9): 2.278342, (11, 9): 2.278342, (12, 9): 0.5695855, (13, 9): 2.278342, (14, 9): 0.14026703, (15, 9): 0.42976263, (11, 10): 2.278342, (12, 10): 0.31445286, (13, 10): 0.14026703, (14, 10): 2.278342, (15, 10): 0.3499277, (12, 11): 0.3162721, (13, 11): 0.42976263, (14, 11): 0.3499277, (15, 11): 2.278342, (13, 12): 2.278342, (14, 12): 2.278342, (15, 12): 2.278342, (14, 13): 2.278342, (15, 13): 2.278342, (15, 14): 2.278342} + + offset: 18.226736272476796 + + + +## Generate a random solution and check its feasibility + + +```python +rng = np.random.default_rng(seed=1234) +random_solution = {i: rng.integers(2) for i in range(F.size)} +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, 8: 0, 9: 0, 10: 1, 11: 0, 12: 1, 13: 0, 14: 1, 15: 0} + + + + +```python +feasibility = qubo_qap_feasibility((F, D), random_solution) +print(f'The feasibility of the random solution is {feasibility}\n') +``` + + The feasibility of the random solution is False + + + +## Generate a feasible solution and check its feasibility + + +```python +feasible_solution = np.zeros(F.shape) +sequence = np.arange(F.shape[0]) +np.random.shuffle(sequence) +for i in range(F.shape[0]): + feasible_solution[i, sequence[i]] = 1 +feasible_solution = {k:v for k, v in enumerate(feasible_solution.flatten())} +print(f'The feasible solution is {feasible_solution}\n') +``` + + The feasible solution is {0: 0.0, 1: 0.0, 2: 1.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 1.0, 8: 0.0, 9: 1.0, 10: 0.0, 11: 0.0, 12: 1.0, 13: 0.0, 14: 0.0, 15: 0.0} + + + + +```python +feasibility = qubo_qap_feasibility((F, D), feasible_solution) +print(f'The feasibility of the feasible solution is {feasibility}\n') +``` + + The feasibility of the feasible solution is True + + + +## Calculate the energy of the feasible solution + + +```python +energy = qubo_qap_energy((F,D), feasible_solution) +print(f'The energy of the feasible solution is {energy}') +``` + + The energy of the feasible solution is 2.7219091992575177 + + +## Hamiltonian + + +```python +ham = hamiltonian_qap((F, D), dense=False) +``` + + [Qibo 0.1.6|INFO|2022-05-31 14:47:26]: Using qibojit backend on /GPU:0 + + +## Solve the Hamiltonian with QAOA + +QAP of size 4 is too large for Qibo QAOA. Let's reduce the size to 3 + + +```python +ham = hamiltonian_qap((F[:3,:3], D[:3,:3]), dense=False) + + +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. diff --git a/examples/qap/main.py b/examples/qap/main.py new file mode 100755 index 0000000000..66146ddf11 --- /dev/null +++ b/examples/qap/main.py @@ -0,0 +1,79 @@ +"""Quadratic Assignment Problem""" + +import argparse + +import numpy as np +from qap import ( + hamiltonian_qap, + qubo_qap, + qubo_qap_energy, + qubo_qap_feasibility, + qubo_qap_penalty, +) +from qubo_utils import binary2spin, spin2QiboHamiltonian + +parser = argparse.ArgumentParser() +parser.add_argument("--filename", default="./tiny04a.dat", type=str) + + +def load_qap(filename): + """Load qap problem from a file + + The file format is compatible with the one used in QAPLIB + + """ + + with open(filename) as fh: + n = int(fh.readline()) + + numbers = [float(n) for n in fh.read().split()] + + data = np.asarray(numbers).reshape(2, n, n) + f = data[1] + d = data[0] + + i = range(len(f)) + f[i, i] = 0 + d[i, i] = 0 + + return f, d + + +def main(filename: str = "./tiny04a.dat"): + print(f"Load flow and distance matrices from {filename} and make a QUBO") + F, D = load_qap(filename) + penalty = qubo_qap_penalty((F, D)) + + linear, quadratic, offset = qubo_qap((F, D), 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 range(F.size)} + feasibility = qubo_qap_feasibility((F, D), random_solution) + assert not feasibility, "The random solution should be infeasible." + + print("Generate a feasible solution and check its feasibility") + feasible_solution = np.zeros(F.shape) + sequence = np.arange(F.shape[0]) + np.random.shuffle(sequence) + for i in range(F.shape[0]): + feasible_solution[i, sequence[i]] = 1 + feasible_solution = {k: v for k, v in enumerate(feasible_solution.flatten())} + feasibility = qubo_qap_feasibility((F, D), feasible_solution) + assert feasibility, "The fixed solution should be feasible." + + print("Calculate the energy of the solution") + energy = qubo_qap_energy((F, D), feasible_solution) + + print("Construct a hamiltonian directly from flow and distance matrices") + ham = hamiltonian_qap((F, D), dense=False) + + print("done.") + + +if __name__ == "__main__": + # by defualt, test on the mvc.csv in the same directory + args = parser.parse_args() + main(args.filename) diff --git a/examples/qap/qap.py b/examples/qap/qap.py new file mode 100644 index 0000000000..1efecae8e9 --- /dev/null +++ b/examples/qap/qap.py @@ -0,0 +1,138 @@ +from typing import Any, Dict, Tuple + +import numpy as np +from qubo_utils import binary2spin, spin2QiboHamiltonian + + +def qubo_qap(g: Tuple[np.ndarray, np.ndarray], penalty: float = None): + """Given adjacency matrix of flow and distance, formulate the QUBO of + Quadratic Assignment Problem (QAP) + + Args: + g: the adjacency matrices of flow and distance that describe the QAP + penalty: penalty weight for the constraints, if not given, it + is inferred from the adjacency matrices + + Returns: + linear (Dict[Any, float]): linear terms + quadratic (Dict[(Any, Any), float]): quadratic terms + offset (float): the constant term + """ + + flow, distance = g + n = len(flow) + q = np.einsum("ij,kl->ikjl", flow, distance) + + if penalty is None: + penalty = qubo_qap_penalty(g) + + i = range(len(q)) + q[i, :, i, :] += 1 * penalty + q[i, :, i, :] -= 2 * penalty * np.eye(n) + q[:, i, :, i] += 1 * penalty + q[:, i, :, i] -= 2 * penalty * np.eye(n) + + q = q.reshape(n**2, n**2).astype(np.float32) + + offset = penalty * 2 * n + linear = {i: q[i, i] for i in range(q.shape[0])} + quadratic = { + (i, j): q[i, j] for j in range(q.shape[1]) for i in range(q.shape[0]) if i > j + } + + return linear, quadratic, offset + + +def qubo_qap_penalty(g: Tuple[np.ndarray, np.ndarray]): + """Find appropriate penalty weight to ensure feasibility + + Args: + g: the adjacency matrices of flow and distance that describe the QAP + + Returns: + penalty (float): penalty weight for the constraints + """ + F, D = g + q = np.einsum("ij,kl->ikjl", F, D) + return F.shape[0] * np.abs(q).max() + + +def qubo_qap_feasibility(g: Tuple[np.ndarray, np.ndarray], solution: Dict[Any, int]): + """Check if the solution violates the constraints of the problem + + Args: + g: the adjacency matrices of flow and distance that describe the QAP + solution (Dict[Any, int]): the binary solution + + Returns: + feasibility (bool): whether the solution meet the constraints + """ + F, D = g + + solution = [[k, v] for k, v in solution.items()] + solution.sort(key=lambda x: x[0]) + _, vals = zip(*solution) + vals = np.asarray(vals).reshape(F.shape) + + assert np.all( + np.logical_or(vals == 0, vals == 1) + ), "Decision variable must be 0 or 1" + return np.all(vals.sum(axis=0) == 1) and np.all(vals.sum(axis=1) == 1) + + +def qubo_qap_energy(g: Tuple[np.ndarray, np.ndarray], solution: Dict[Any, int]): + """Calculate the energy of the solution on a QAP problem + The result is based on the assumption that soltuion is feasible. + + Args: + g: the adjacency matrices of flow and distance that describe the QAP + solution (Dict[Any, int]): the solution + + Returns: + energy (float): the energy of the solution + """ + + F, D = g + + solution = [[k, v] for k, v in solution.items()] + solution.sort(key=lambda x: x[0]) + _, vals = zip(*solution) + vals = np.asarray(vals).reshape(F.shape) + + state = np.vstack(np.where(vals == 1)).T + state = np.array(state).astype(np.int32) + energy = 0 + for i in range(len(state)): + energy += np.einsum( + "i,i->", F[state[i, 0], state[:, 0]], D[state[i, 1], state[:, 1]] + ) + return energy + + +def hamiltonian_qap( + g: Tuple[np.ndarray, np.ndarray], penalty: float = None, dense: bool = False +): + """Given a flow and distance matrices, return the hamiltonian of Quadratic + Assignment Problem (QAP). + + If penalty is not given, it is inferred from g. + + Args: + g (Tuple[np.ndarray, np.ndarray]): flow and distance matrices + penalty (float): penalty weight of the constraints + dense (bool): sparse or dense hamiltonian + + Returns: + ham: qibo hamiltonian + + """ + + linear, quadratic, offset = qubo_qap(g, penalty) + + 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, dense=dense) + + return ham diff --git a/examples/qap/qubo_utils.py b/examples/qap/qubo_utils.py new file mode 100644 index 0000000000..33abe53057 --- /dev/null +++ b/examples/qap/qubo_utils.py @@ -0,0 +1,122 @@ +from typing import Dict, Tuple + + +def binary2spin( + linear: Dict[int, float], + quadratic: Dict[Tuple[int, int], float], + offset: float = 0, +): + """Convert binary model to spin model + + Please remember to put a negative sign to h and J after using this + function , if you are going to form a mamiltonian from the spin + model. Hamiltonians usually have a leading negative sign in them, + but QUBOs don't. + + Args: + linear (dict): linear term of the binary model + quadratic (dict): quadratic term of the binary model + offset (float): offset of the binary model + + Returns: + h (dict): bias of the spin model + J (dict): interaction of the spin model + offset (float): offset of the spin model + """ + + h = {x: 0.5 * w for x, w in linear.items()} + + J = [] + for (x, y), w in quadratic.items(): + J.append(((x, y), 0.25 * w)) + h[x] += 0.25 * w + h[y] += 0.25 * w + J = dict(J) + + offset += 0.5 * sum(linear.values()) + offset += 0.25 * sum(quadratic.values()) + + return h, J, offset + + +def spin2binary( + h: Dict[int, float], + J: Dict[Tuple[int, int], float], + offset: float = 0, +): + """Convert spin model to binary model + + Please remember to put a negative sign to h and J before using this + function if you extract them from a hamiltonian. Hamiltonians + usually have a leading negative sign in them, but QUBOs don't. + + Args: + h (dict): bias of the spin model + J (dict): interaction of the spin model + offset (float): offset of the spin model + + Returns: + linear (dict): linear term of the binary model + quadratic (dict): quadratic term of the binary model + offset (float): offset of the binary model + """ + + linear = {s: 2.0 * bias for s, bias in h.items()} + + quadratic = [] + for (s, t), bias in J.items(): + quadratic.append(((s, t), 4.0 * bias)) + linear[s] -= 2.0 * bias + linear[t] -= 2.0 * bias + quadratic = dict(quadratic) + + offset -= sum(linear.values()) + offset += sum(quadratic.values()) + + return linear, quadratic, offset + + +def spin2QiboHamiltonian( + h: Dict[int, float], + J: Dict[Tuple[int, int], float], + dense: bool = True, +): + """Convert spin model to qibo Hamiltonian + + Mixer is not included. + + Please remember to put a negative sign to h and J if you get h and J + from binary2spin. Hamiltonians usually have a leading negative sign + in them, but QUBOs don't. + + .. math:: + H = - \\sum_{i=0}^N \\sum _{j=0}^N J_{ij} Z_i Z_j - \\sum_{i=0}^N h_i Z_i. + + Args: + h (dict): bias of the spin model + J (dict): interaction of the spin model + dense (bool): If ``True`` it creates the Hamiltonian as a + :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates + a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`. + + Returns: + linear (dict): linear term of the binary model + quadratic (dict): quadratic term of the binary model + offset (float): offset of the binary model + """ + + from qibo import hamiltonians + from qibo.symbols import Z + + form = 0 + for k, v in h.items(): + form -= Z(k, commutative=True) * v + for (k0, k1), v in J.items(): + form -= Z(k0, commutative=True) * Z(k1, commutative=True) * v + + ham = hamiltonians.SymbolicHamiltonian(form) + + if dense: + return ham.dense + else: + return ham diff --git a/examples/qap/tiny04a.dat b/examples/qap/tiny04a.dat new file mode 100644 index 0000000000..39661414a2 --- /dev/null +++ b/examples/qap/tiny04a.dat @@ -0,0 +1,11 @@ +4 + +0.0 0.7796977849571155 0.4304502176328382 0.4329405547124155 +0.7796977849571155 0.0 0.1920096021905386 0.5882961822495594 +0.4304502176328382 0.1920096021905386 0.0 0.47901121613954656 +0.4329405547124155 0.5882961822495594 0.47901121613954656 0.0 + +0.0 0.29541330831110774 0.6844285463729066 0.1988227908746449 +0.29541330831110774 0.0 0.6164922460559509 0.16210678581988341 +0.6844285463729066 0.6164922460559509 0.0 0.7305208755290076 +0.1988227908746449 0.16210678581988341 0.7305208755290076 0.0