From dba4072e6d46070e44965a3310add539a938181c Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Thu, 13 Jul 2023 16:55:35 -0700 Subject: [PATCH] Add first VQE submodule iteration Co-authored-by: Juan Sebastian Lozano --- dwave/gate/operations/operations.py | 21 + dwave/gate/vqe/__init__.py | 21 + dwave/gate/vqe/cqm_partition.py | 68 ++ dwave/gate/vqe/operators.py | 624 ++++++++++++++++++ .../gate/vqe/transformations/jordan_wigner.py | 170 +++++ dwave/gate/vqe/vqe.py | 136 ++++ 6 files changed, 1040 insertions(+) create mode 100644 dwave/gate/vqe/__init__.py create mode 100644 dwave/gate/vqe/cqm_partition.py create mode 100644 dwave/gate/vqe/operators.py create mode 100644 dwave/gate/vqe/transformations/jordan_wigner.py create mode 100644 dwave/gate/vqe/vqe.py diff --git a/dwave/gate/operations/operations.py b/dwave/gate/operations/operations.py index 93cf509..05a85e7 100644 --- a/dwave/gate/operations/operations.py +++ b/dwave/gate/operations/operations.py @@ -29,6 +29,7 @@ "Y", "Z", "Hadamard", + "HY", "Phase", "S", # alias "P", # alias @@ -231,6 +232,26 @@ def matrix(cls) -> NDArray[np.complex128]: matrix = math.sqrt(2) / 2 * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex128) return matrix +class HY(Operation): + """Self-inverse operation for Y-Z basis swap. + + Has the property $\sigma_y = H_y \sigma_z H_y. + + Args: + qubits: Qubits on which the operation should be applied. Only + required when applying an operation within a circuit context. + """ + + _num_qubits: int = 1 + + # see: https://quantumcomputing.stackexchange.com/questions/11390/problem-with-commutation-of-e-ih-1t-and-e-ih-2t-where-h-1-commutes/11391#11391 + # not the only way of doing this, but the most simple + + @mixedproperty + def matrix(cls) -> NDArray[np.complex128]: + """The matrix representation of the $H_y$ operator.""" + return math.sqrt(2) / 2 * np.array([[1, -1j], [1j, -1]], dtype=np.complex128) + class S(Operation): """S (Phase) operation. diff --git a/dwave/gate/vqe/__init__.py b/dwave/gate/vqe/__init__.py new file mode 100644 index 0000000..8c836c2 --- /dev/null +++ b/dwave/gate/vqe/__init__.py @@ -0,0 +1,21 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Module for running the Variational Quantum Eigensolver (VQE) algorithm. + +TODO: Add description +""" +from dwave.gate.vqe.operators import * +from dwave.gate.vqe.transformations.jordan_wigner import * +from dwave.gate.vqe.vqe import * diff --git a/dwave/gate/vqe/cqm_partition.py b/dwave/gate/vqe/cqm_partition.py new file mode 100644 index 0000000..cbfe9cf --- /dev/null +++ b/dwave/gate/vqe/cqm_partition.py @@ -0,0 +1,68 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +__all__ = [ + "partition_cqm", +] + +import warnings + + +try: + import dimod + + from dwave.system import LeapHybridCQMSampler +except ImportError as e: + raise ImportError("Dimod must be installed for CQM partioning.") from e + + +def partition_cqm(adjacency, num_measurements): + """TODO""" + cqm = dimod.ConstrainedQuadraticModel() + + measurement_rounds = [ + dimod.Integer(i, lower_bound=0, upper_bound=num_measurements) + for i in range(num_measurements) + ] + + cqm.set_objective(sum(measurement_rounds)) + + for source, target in adjacency: + cqm.add_constraint_from_comparison( + measurement_rounds[source] ** 2 + + measurement_rounds[target] ** 2 + - 2 * measurement_rounds[source] * measurement_rounds[target] + >= 1 + ) + + sampler = LeapHybridCQMSampler() + + sampleset = sampler.sample_cqm(cqm) + sampleset = sampleset.filter(lambda x: x.is_feasible) + + if len(sampleset) == 0: + warnings.warn("No feasible measurement schemes found, defaulting to sequential measurement") + cqm_sample = {i: i for i in range(num_measurements)} + else: + cqm_sample = sampleset.first.sample + + partition = {} + + for pauli_index, measurement_round in cqm_sample.items(): + if measurement_round not in partition: + partition[measurement_round] = [pauli_index] + else: + partition[measurement_round] = partition[measurement_round] + [pauli_index] + + return partition diff --git a/dwave/gate/vqe/operators.py b/dwave/gate/vqe/operators.py new file mode 100644 index 0000000..8950ee8 --- /dev/null +++ b/dwave/gate/vqe/operators.py @@ -0,0 +1,624 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from __future__ import annotations + +__all__ = [ + "tensor", + "Tensor", + "Pauli", + "OperatorTerm", + "Operator", +] + +import copy +import numbers +from itertools import product +from typing import Any, List, Sequence, Tuple, Union + +import numpy as np + +def tensor(*factors: Sequence[Pauli]) -> Union[Tensor, Pauli]: + """TODO""" + if len(factors) == 0: + return Pauli() + if len(factors) == 1: + return factors[0] + + return Tensor(*factors) + + +class Tensor: + """TODO""" + + def __init__(self, *factors: Sequence[Pauli]) -> None: + if len(factors) < 2: + raise ValueError("Minimum of 2 Pauli objects required to create a tensor.") + for f in factors: + if not isinstance(f, Pauli): + raise TypeError("Only Pauli object are accepted factors.") + + # store as tuple to make immutable/hashable + self._factors = tuple(factors) + + @property + def indices(self) -> List[int]: + """TODO""" + return [f.idx for f in self.factors] + + @property + def factors(self) -> Tuple[Pauli, ...]: + """TODO""" + return self._factors + + @property + def is_identity(self) -> bool: + """TODO""" + return all(f.is_identity for f in self.factors) + + def __matmul___(self, pauli: Tensor) -> Tensor: + return Tensor(*self.factors, pauli) + + def __len__(self) -> int: + return len(self.factors) + + def __getitem__(self, item: int) -> Pauli: + return self.factors[item] + + def __reversed__(self) -> Union[Tensor, Pauli]: + if isinstance(self, Pauli): + return self + return Tensor(*self.factors[::-1]) + + def __eq__(self, tensor: Tensor) -> bool: + if isinstance(tensor, Tensor) and self.factors == tensor.factors: + return True + return False + + def __hash__(self) -> int: + return hash(self.factors) + + def __repr__(self) -> str: + return str(self) + + def __str__(self) -> str: + return " @ ".join(map(str, self.factors)) + + +class Pauli(Tensor): + """TODO""" + + X = np.array([[0, 1], [1, 0]], dtype=complex) + Y = np.array([[0, complex(0, 1)], [-complex(0, 1), 0]], dtype=complex) + Z = np.array([[1, 0], [0, -1]], dtype=complex) + Id = np.array([[1, 0], [0, 1]], dtype=complex) + + pauli_algebra = {"X": X, "Y": Y, "Z": Z, "I": Id} + + def __init__(self, name=None, idx=None) -> None: + if name and name not in Pauli.pauli_algebra.keys(): + raise ValueError( + f"Name not in pauli algebra, must be one of {Pauli.pauli_algebra.keys()}" + ) + + self.name = name or "I" + self.idx = idx + self.value = Pauli.pauli_algebra[self.name] + self._factors = [self] + + @property + def is_identity(self) -> bool: + """TODO""" + return self.name == "I" + + @staticmethod + def match_pauli_matrix(p: Pauli) -> Tuple[str, complex, List]: + """TODO""" + if np.equal(np.identity(2, dtype=complex), p).all(): + return 1 + + matches = [key for key, val in Pauli.pauli_algebra.items() if np.equal(val, p).all()] + + if len(matches) == 1: + return matches[0] + + matches = [ + key + for key, val in Pauli.pauli_algebra.items() + if np.equal(val, p * complex(0, 1)).all() + ] + if len(matches) == 1: + return ("*", complex(0, 1), matches[0]) + + matches = [ + key for key, val in Pauli.pauli_algebra.items() if np.equal(val, -1 * p).all() + ] + if len(matches) == 1: + return ("*", -1, matches[0]) + + matches = [ + key + for key, val in Pauli.pauli_algebra.items() + if np.equal(val, p * complex(0, -1)).all() + ] + if len(matches) == 1: + return ("*", complex(0, -1), matches[0]) + + + def __matmul___(self, pauli: Pauli) -> Tensor: + return Tensor(self, pauli) + + def __mul__(self, other: Any) -> Tuple[str, Pauli, Any]: + if isinstance(other, Pauli): + return Pauli.match_pauli_matrix(self.value @ other.value) + else: + return ("*", self, other) + + def __eq__(self, pauli: Pauli) -> bool: + if isinstance(pauli, Pauli) and self.name == pauli.name and self.idx == pauli.idx: + return True + return False + + def __hash__(self) -> int: + return hash(self.name) + hash(self.idx) + + def __str__(self) -> str: + if self.idx is not None: + return f"{self.name}{self.idx}" + return f"{self.name}" + + def __repr__(self) -> str: + return str(self) + + +class OperatorTerm: + """TODO""" + + supported_operators = ["*", "+"] + + def __init__(self, operator: str, left, right) -> None: + if operator not in OperatorTerm.supported_operators: + raise ValueError( + f"Operator not supported, must be one of {OperatorTerm.supported_operators}" + ) + + self.operator = operator + self.left = left + self.right = right + + def reduce(self): + """TODO""" + if isinstance(self.left, OperatorTerm): + self.left = self.left.reduce() + + if isinstance(self.right, OperatorTerm): + self.right = self.right.reduce() + + # distributive law + if self.operator == "*": + # a*(x + y) => a*x + a*y + if isinstance(self.right, OperatorTerm) and self.right.operator == "+": + self.operator = "+" + child_left = self.left + self.left = OperatorTerm("*", child_left, self.right.left) + self.right = OperatorTerm("*", child_left, self.right.right) + # (a + b)*x => a*x + b*x + elif isinstance(self.left, OperatorTerm) and self.left.operator == "+": + self.operator = "+" + child_right = self.right + # must be done in this order to not mess up the state + self.right = OperatorTerm("*", self.left.right, child_right) + self.left = OperatorTerm("*", self.left.left, child_right) + + # simplification of the same type + if isinstance(self.left, numbers.Number) and isinstance(self.right, numbers.Number): + if self.operator == "*": + return self.left * self.right + if self.operator == "+": + return self.left + self.right + elif isinstance(self.left, np.ndarray) and isinstance(self.right, np.ndarray): + if self.operator == "*": + return self.left @ self.right + if self.operator == "+": + return self.left + self.right + elif isinstance(self.left, Pauli) and isinstance(self.right, Pauli): + if self.operator == "*": + return_value = self.left * self.right + if isinstance(return_value, numbers.Number): + return return_value + else: + return OperatorTerm(*return_value) + + # identity simplification + if self.operator == "*": + if isinstance(self.left, numbers.Number) and self.left == 0: + return 0 + elif isinstance(self.right, numbers.Number) and self.right == 0: + return 0 + if isinstance(self.left, numbers.Number) and self.left == 1: + return self.right + elif isinstance(self.right, numbers.Number) and self.right == 1: + return self.left + elif self.operator == "+": + if isinstance(self.left, numbers.Number) and self.left == 0: + return self.right + elif isinstance(self.right, numbers.Number) and self.right == 0: + return self.left + + return self + + def split(self): + """TODO""" + split_terms = [] + + if self.operator == "+": + if isinstance(self.left, OperatorTerm): + split_terms.extend(self.left.split()) + else: + split_terms.append(self.left) + + if isinstance(self.right, OperatorTerm): + split_terms.extend(self.right.split()) + else: + split_terms.append(self.right) + + else: + split_terms = [self] + + return split_terms + + def iter_left(self, order="bottom->top"): + """TODO""" + orders = ["bottom->top", "top->bottom"] + + if order not in orders: + raise ValueError(f"method must be one of {orders}, got {order}") + + if order == "top->bottom": + yield self.left + + if isinstance(self.left, OperatorTerm): + yield from self.left.iter_left() + + if order == "bottom->top": + yield self.left + + def iter_right(self, order="bottom->top"): + """TODO""" + orders = ["bottom->top", "top->bottom"] + + if order not in orders: + raise ValueError(f"method must be one of {orders}, got {order}") + + if order == "top->bottom": + yield self.right + + if isinstance(self.left, OperatorTerm): + yield from self.right.iter_right() + + if order == "bottom->top": + yield self.right + + def traverse(self, order="left->right"): + """TODO""" + orders = ["left->right", "right->left"] + + if order not in orders: + raise ValueError(f"method must be one of {orders}, got {order}") + + yield self + + if order == "left->right": + if isinstance(self.left, OperatorTerm): + yield from self.left.traverse(order) + + if isinstance(self.right, OperatorTerm): + yield from self.right.traverse(order) + + elif order == "right->left": + if isinstance(self.right, OperatorTerm): + yield from self.right.traverse(order) + + if isinstance(self.left, OperatorTerm): + yield from self.left.traverse(order) + + def iter_leaves(self, order="left->right"): + """TODO""" + orders = ["left->right", "right->left"] + + if order not in orders: + raise ValueError(f"method must be one of {orders}, got {order}") + + if order == "left->right": + if isinstance(self.left, OperatorTerm): + yield from self.left.iter_leaves(order) + else: + yield self.left + + if isinstance(self.right, OperatorTerm): + yield from self.right.iter_leaves(order) + else: + yield self.right + + elif order == "right->left": + if isinstance(self.right, OperatorTerm): + yield from self.right.iter_leaves(order) + else: + yield self.right + + if isinstance(self.left, OperatorTerm): + yield from self.left.iter_leaves(order) + else: + yield self.left + + def homogenize(self, index=0): + """TODO""" + copy_self = copy.copy(self) + + substituion_dict = {} + substituion_dict_left = {} + substituion_dict_right = {} + + if isinstance(self.left, OperatorTerm): + if self.left.operator == self.operator: + copy_self.left, substituion_dict_left = self.left.homogenize( + index + 1 + ) # indexes arent unique, a bug! + else: + substituion_dict[f"X_{index}"] = self.left + copy_self.left = f"X_{index}" + index += 1 + + if isinstance(self.right, OperatorTerm): + if self.right.operator == self.operator: + copy_self.right, substituion_dict_right = self.right.homogenize(index + 1) + else: + substituion_dict[f"X_{index}"] = self.right + copy_self.right = f"X_{index}" + index += 1 + + # NOTE: Only Python >=3.9 + # return copy_self, substituion_dict | substituion_dict_left | substituion_dict_right + return copy_self, {**substituion_dict, **substituion_dict_left, **substituion_dict_right} + + def collect(self): + """TODO""" + if self.operator == "+": + unit = 0 + elif self.operator == "*": + unit = 1 + + old_term_homogenized, substitutions = self.homogenize() + + final_terms = {} + for term in old_term_homogenized.iter_leaves(): + if isinstance(term, numbers.Number): + if "scalar" not in final_terms: + final_terms["scalar"] = unit + + if self.operator == "+": + final_terms["scalar"] = final_terms["scalar"] + term + elif self.operator == "*": + final_terms["scalar"] = final_terms["scalar"] * term + + else: + if str(term.__class__) not in final_terms: + if term in substitutions: + term = substitutions[term].collect() + final_terms[str(term.__class__)] = term + elif isinstance(term, np.ndarray): + # there is a slight bug here + # I'm assuming that the numpy matrix commutes with the Pauli matracies + # which in general isn't true + if self.operator == "+": + final_terms[str(term.__class__)] = final_terms[str(term.__class__)] + term + elif self.operator == "*": + final_terms[str(term.__class__)] = final_terms[str(term.__class__)] * term + elif isinstance(term, Pauli): + if self.operator == "+": + final_terms[str(term.__class__)] = OperatorTerm( + "+", final_terms[str(term.__class__)], term + ) + elif self.operator == "*": + prod = final_terms[str(term.__class__)] * term + if prod == 1: + del final_terms[str(term.__class__)] + elif len(prod) > 1 and prod[0] == "*": + # this section has some assumptions about the output type of PauliMatrix multiplication + if isinstance(prod[1], numbers.Number): + final_terms["scalar"] = final_terms["scalar"] * prod[1] + elif isinstance(prod[1], Pauli): + final_terms[str(term.__class__)] = prod[1] + + if isinstance(prod[2], numbers.Number): + final_terms["scalar"] = final_terms["scalar"] * prod[2] + elif isinstance(prod[2], Pauli): + final_terms[str(term.__class__)] = prod[2] + else: + if term in substitutions: + term = substitutions[term].collect() + final_terms[str(term.__class__)] = OperatorTerm( + self.operator, final_terms[str(term.__class__)], term + ) + + final_item = unit + + for _, item in final_terms.items(): + final_item = OperatorTerm(self.operator, final_item, item) + + return final_item + + def __mul__(self, other): + return OperatorTerm("*", self, other) + + def __add__(self, other): + return OperatorTerm("+", self, other) + + def __str__(self) -> str: + return "(" + str(self.operator).join([str(self.left), str(self.right)]) + ")" + + def __repr__(self) -> str: + return str(self) + + def __eq__(self, other): + # super strict version of equality - for example non-associative + # subsuming some errors in favor of easy typing + try: + return (self.left == other.left) and (self.right == other.right) + except Exception as e: + return False + + +class Operator: + """TODO""" + + def __init__(self, length: int = None, values=None, default_value=0) -> None: + if (values is None) and (length is None): + raise ValueError(f"Either length of values must be specified") + + if values is None: + self.operator_terms = [default_value for _ in range(length)] + elif (length is None) or (len(values) == length): + self.operator_terms = values + else: + raise ValueError(f"Wrong length of values: given {len(values)}, expected {length}") + + self.length = len(self.operator_terms) + + def reduce(self): + """TODO""" + for index, term in enumerate(self.operator_terms): + if isinstance(term, OperatorTerm): + reduced = False + while not reduced: + old_term = copy.copy(self.operator_terms[index]) + self.operator_terms[index] = self.operator_terms[index].reduce() + if ( + not isinstance(self.operator_terms[index], OperatorTerm) + or old_term == self.operator_terms[index] + ): + reduced = True + if isinstance(self.operator_terms[index], OperatorTerm): + self.operator_terms[index] = self.operator_terms[index].collect() + + if isinstance(self.operator_terms[index], OperatorTerm): + self.operator_terms[index] = self.operator_terms[index].reduce() + + if not reduced and not isinstance(self.operator_terms[index], OperatorTerm): + reduced = True + + def split(self): + """TODO""" + return_values = [] + for index in range(self.length): + if isinstance(self.operator_terms[index], OperatorTerm): + return_values.append(self.operator_terms[index].split()) + else: + return_values.append([self.operator_terms[index]]) + + return [Operator(values=i) for i in product(*return_values)] + + def collect(self): + """TODO""" + for index in range(self.length): + if isinstance(self.operator_terms[index], OperatorTerm): + self.operator_terms[index] = self.operator_terms[index].collect() + + def to_pauli_strings(self): + """TODO""" + split_operators = self.split() + + pauli_strings = [] + for operator in split_operators: + coefficient = 1 + paulis = [] + for index, term in enumerate(operator.operator_terms): + if isinstance(term, numbers.Number): + coefficient *= term + elif isinstance(term, Pauli): + paulis.append(Pauli(term.name, index)) + elif isinstance(term, OperatorTerm): + if not ( + isinstance(term.left, numbers.Number) or isinstance(term.left, Pauli) + ): + raise ValueError( + f"term left child is of unsupported class {term.left.__class__}, likely not simplified enough" + ) + + if not ( + isinstance(term.right, numbers.Number) + or isinstance(term.right, Pauli) + ): + raise ValueError( + f"term right child is of unsupported class {term.left.__class__}, likely not simplified enough" + ) + + for child in [term.left, term.right]: + if isinstance(child, numbers.Number): + coefficient *= child + elif isinstance(child, Pauli): + paulis.append(Pauli(child.name, index)) + else: + raise ValueError( + f"term {term} of class {term.__class__} not a Number, PauliMatrix, or OperatorTerm" + ) + pauli_strings.append((coefficient, tensor(*paulis))) + return pauli_strings + + def __mul__(self, other): + if not isinstance(other, Operator): + raise ValueError( + f"Operators can only be multiplied with operators: got {other.__class__}" + ) + + if self.length != other.length: + raise ValueError( + f"Operators must be of same length: got {self.length} and {other.length}" + ) + + return Operator( + values=[ + OperatorTerm("*", i, j) for i, j in zip(self.operator_terms, other.operator_terms) + ] + ) + + def _add__(self, other): + if not isinstance(other, Operator): + raise ValueError( + f"Operators can only be multiplied with operators: got {other.__class__}" + ) + + if self.length != other.length: + raise ValueError( + f"Operators must be of same length: got {self.length} and {other.length}" + ) + + return Operator( + values=[ + OperatorTerm("+", i, j) for i, j in zip(self.operator_terms, other.operator_terms) + ] + ) + + def __str__(self) -> str: + return "$" + " \otimes ".join([str(i) for i in self.operator_terms]) + "$" + + def __repr__(self) -> str: + return str(self) + + def __eq__(self, other): + if not isinstance(other, Operator): + raise ValueError( + f"Operators can only be multiplied with operators: got {other.__class__}" + ) + + return all([i == j for i, j in zip(self.operator_terms, other.operator_terms)]) diff --git a/dwave/gate/vqe/transformations/jordan_wigner.py b/dwave/gate/vqe/transformations/jordan_wigner.py new file mode 100644 index 0000000..094da1a --- /dev/null +++ b/dwave/gate/vqe/transformations/jordan_wigner.py @@ -0,0 +1,170 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +__all__ = [ + "JordanWignerAnnihilation", + "JordanWignerCreation", + "JordanWignerHamiltonian", +] + +from itertools import product + +import numpy as np + +from dwave.gate.vqe.operators import Operator, OperatorTerm, Pauli + + +class JordanWignerAnnihilation(Operator): + """TODO""" + + def __init__(self, index: int, length: int = None) -> None: + values = [] + for _ in range(index): + values.append(Pauli("Z")) + + values.append( + OperatorTerm( + "*", + 0.5, + OperatorTerm( + "+", Pauli("X"), OperatorTerm("*", complex(0, 1), Pauli("Y")) + ), + ) + ) + + for _ in range(index, length): + values.append(1) + + super().__init__(values=values) + + +class JordanWignerCreation(Operator): + """TODO""" + + def __init__(self, index: int, length: int = None) -> None: + values = [] + for _ in range(index): + values.append(Pauli("Z")) + + values.append( + OperatorTerm( + "*", + 0.5, + OperatorTerm( + "+", Pauli("X"), OperatorTerm("*", complex(0, -1), Pauli("Y")) + ), + ) + ) + + for _ in range(index, length): + values.append(1) + + super().__init__(values=values) + + +class JordanWignerHamiltonian: + """TODO""" + + def __init__( + self, num_orbitals: int, first_order: np.ndarray, second_order: np.ndarray + ) -> None: + if (len(first_order.shape) != 2) or any([dim != num_orbitals for dim in first_order.shape]): + raise ValueError( + "First order interactions of wrong shape, expected shape " + f"{(num_orbitals, num_orbitals)}, got {first_order.shape}" + ) + + if (len(second_order.shape) != 4) or any( + [dim != num_orbitals for dim in second_order.shape] + ): + raise ValueError( + "Second order interactions of wrong shape, expected shape " + f"{(num_orbitals, num_orbitals, num_orbitals, num_orbitals)}, " + f"got {second_order.shape}" + ) + + self.num_orbitals = num_orbitals + + self.terms = {} + + for idx_1, idx_2 in product(range(num_orbitals), repeat=2): + a = JordanWignerAnnihilation(idx_1, num_orbitals) + a_dagger = JordanWignerCreation(idx_2, num_orbitals) + + term = a_dagger * a + + term.reduce() + + self.terms[(idx_1, idx_2)] = (first_order[idx_1, idx_2], term.to_pauli_strings()) + + for idx_1, idx_2, idx_3, idx_4 in product(range(num_orbitals), repeat=4): + a_1 = JordanWignerAnnihilation(idx_1, num_orbitals) + a_2 = JordanWignerAnnihilation(idx_2, num_orbitals) + a_1_dagger = JordanWignerCreation(idx_3, num_orbitals) + a_2_dagger = JordanWignerCreation(idx_4, num_orbitals) + + term = a_1_dagger * a_1 * a_2 * a_2_dagger + + term.reduce() + + self.terms[(idx_1, idx_2, idx_3, idx_4)] = ( + second_order[idx_1, idx_2, idx_3, idx_4], + term.to_pauli_strings(), + ) + + self.pauli_terms = None + + def collect_pauli_terms(self): + """TODO""" + if self.pauli_terms is None: + # assumes commutivity of sum terms (not sure if okay?) + sum_of_coef = {} + for item in self.terms.values(): + classical_coef = item[0] + for qubit_coef, pauli_string in item[1]: + + if pauli_string.is_identity: + continue + + if pauli_string not in sum_of_coef: + sum_of_coef[pauli_string] = 0 + + sum_of_coef[pauli_string] += classical_coef * qubit_coef + self.pauli_terms = sum_of_coef + + return self.pauli_terms + + def calculate_energy(self, measurements): + """TODO""" + energy = 0 + for item in self.terms.values(): + classical_coef = item[0] + for qubit_coef, pauli_string in item[1]: + if not pauli_string.is_identity and (pauli_string not in measurements): + raise ValueError(f"No measurement found for {pauli_string}") + elif pauli_string.is_identity: + energy += classical_coef * qubit_coef + else: + energy += classical_coef * qubit_coef * measurements[pauli_string] + + return energy + + def __repr__(self) -> str: + return ( + f"Jordan-Wigner molecular hamiltonian with {self.num_orbitals} " + f"orbitals and {len(self.terms)} hamiltonian terms" + ) + + def __str__(self) -> str: + return self.__repr__() diff --git a/dwave/gate/vqe/vqe.py b/dwave/gate/vqe/vqe.py new file mode 100644 index 0000000..31945c1 --- /dev/null +++ b/dwave/gate/vqe/vqe.py @@ -0,0 +1,136 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +__all__ = [ + "initial_parameters", + "Measurements", +] + +import warnings + +import numpy as np + +from dwave.gate.vqe.transformations.jordan_wigner import JordanWignerHamiltonian + +try: + from dwave.gate.vqe.cqm_partition import partition_cqm +except ImportError: + cqm_supported = False +else: + from dwave.cloud.exceptions import SolverFailureError + + cqm_supported = True + + +def initial_parameters(hamiltonian: JordanWignerHamiltonian): + """TODO""" + pauli_terms = hamiltonian.collect_pauli_terms() + + for string, coef in pauli_terms.items(): + if coef.imag != 0: + # I have no idea what to do with the imaginary coefficients + # For the trotterization we need only real coefficients + # For now I take the norm, but this means something is either (1) + # broken or (2) the ansatz differs from the literature in Anand et + # al 2022 they use a slightly different ansatz but they say there + # shouldn't be imag coefficients..... + pauli_terms[string] = np.exp(coef).real + else: + pauli_terms[string] = np.exp(coef) + + # no idea if this is right + return [val for _, val in pauli_terms.items()] + + +class Measurements: + """TODO""" + + def __init__(self, hamiltonian: JordanWignerHamiltonian) -> None: + self.pauli_terms = list(hamiltonian.collect_pauli_terms().keys()) + + self.num_qubits = hamiltonian.num_orbitals + self.num_measurements = len(self.pauli_terms) + self.adjacency = None + + def generate_adjacency(self): + """TODO""" + self.adjacency = [] + for index, term in enumerate(self.pauli_terms): + for previous_index, previous_term in enumerate(self.pauli_terms[:index]): + commute_sum = self.num_qubits - len(term) # commute on the identity + for t in term.factors: + if (t.idx not in previous_term.indices) or previous_term[t.idx] == t: + commute_sum += 1 + if commute_sum % 2 == 0: + # two Pauli strings do not commute if and only if they do + # commute on an even number of indices + self.adjacency.append((index, previous_index)) + + def _partition_linear(self): + """TODO""" + measurement_assignments = {} + + current_measurement = 0 + for term in range(self.num_measurements): + new_measurement = any( + [ + ((i, term) in self.adjacency) or ((term, i) in self.adjacency) + for i in range(term) + ] + ) + if new_measurement: + current_measurement += 1 + + measurement_assignments[term] = current_measurement + + partition = {} + + for pauli_index, measurement_round in measurement_assignments.items(): + if measurement_round not in partition: + partition[measurement_round] = [pauli_index] + else: + partition[measurement_round] = partition[measurement_round] + [pauli_index] + + return partition + + def partition(self, method: str = "best"): + """TODO""" + supported_methods = ["best", "cqm", "linear"] + + if method not in supported_methods: + raise ValueError(f"method must be one of {supported_methods}, got {method}") + + if self.adjacency is None: + self.generate_adjacency() + + if method == "cqm" or (method == "best" and cqm_supported): + if not cqm_supported: + raise ValueError("Cannot use CQM method unless dimod is installed.") + + try: + partition = partition_cqm(self.adjacency, self.num_measurements) + except SolverFailureError as e: + # only proceed with 'linear' if 'method' is 'best' + if method != "best": + raise e + + warnings.warn(str(e) + " Using 'linear' method instead.", stacklevel=2) + method = "linear" + else: + return partition + + if method == "linear" or (method == "best" and not cqm_supported): + partition = self._partition_linear() + + return partition