Add the ability to load Pauli strings

qre/qre_hamiltonian.py

@@ -27,16 +27,68 @@ def boson_to_qubit_operator(bosonic_operator, Nmax):
 
 # -------------------------------------------------------------------------------------------------
 
-# TODO: This is a generally useful utility.  Where should it live?
-def tuple_to_string(pauli_tuple, coef, num_qubits):
-    s = ["I",] * num_qubits
-    for idx, op in pauli_tuple:
-        s[idx] = op
-    s = "".join(s)
-    return (s, coef)
+# TODO: These are generally useful utilities.  Where should they live?
+def sparse_to_dense_pauli(sparse_pauli, num_qubits):
+    dense_pauli = ["I",] * num_qubits
+    for idx, op in sparse_pauli:
+        dense_pauli[idx] = op
+    return "".join(dense_pauli)
+def dense_to_sparse_pauli(dense_pauli):
+    sparse_pauli = tuple()
+    for idx, op in enumerate(dense_pauli):
+        if op in ["X", "Y", "Z"]:
+            sparse_pauli = (*sparse_pauli, (idx,op))
+        elif op != "I":
+            raise ValueError(f"Invalid character in dense pauli string: \"{op}\".")
+    return sparse_pauli
 
 # -------------------------------------------------------------------------------------------------
 
+class LinearCombinationOfPauliStrings:
+    def __init__(self, **kwargs):
+        self._nq = None
+        self._format = None
+        self._data = None
+        self._nq = kwargs["num_qubits"]
+        for f in [ "dense", "sparse" ]:
+            if f in kwargs:
+                if self._format is not None:
+                    raise ValueError(
+                        "Too many formats provided to LinearCombinationOfPauliStrings.")
+                self._format = f
+                self._data = kwargs[f]
+                assert isinstance(self._data, dict)
+        if self._format is None:
+            raise ValueError("No data provided to LinearCombinationOfPauliStrings.")
+    def num_qubits(self):
+        raise NotImplementedError()
+    def get_dense_pauli_strings(self):
+        if self._format == "dense":
+            return self._data
+        elif self._format == "sparse":
+            return {sparse_to_dense_pauli(pauli, self._nq) : coef
+                    for pauli, coef in self._data.items()}
+        else:
+            raise ValueError("Invalid data format \"{self._format}\".")
+    def get_sparse_pauli_strings(self):
+        if self._format == "dense":
+            return {dense_to_sparse_pauli(pauli) : coef for pauli, coef in self._data.items()}
+        elif self._format == "sparse":
+            return self._data
+        else:
+            raise ValueError("Invalid data format \"{self._format}\".")
+    def energy_shift(self, shift):
+        all_identity = tuple()
+        if self._format == "dense":
+            all_identity = sparse_to_dense_pauli(all_identity, self._nq)
+        identity_coefficient = self._data.get(all_identity, 0.0) + shift
+        self._data[all_identity] = identity_coefficient
+
+# -------------------------------------------------------------------------------------------------
+
+# TODO: The heavy use of isinstance() suggests that perhaps Hamiltonian should be a base class that
+#       other things are built on top of?
+
 class Hamiltonian:
     def __init__(self, hamiltonian):
         self._H = hamiltonian
@@ -45,13 +97,15 @@ def get_core_operator(self):
         #       pyLIQTR problem instance.
         return self._H
     def set_fermionic_mapping(self, mapping):
+        # self._fmap is never used for LinearCombinationOfPauliStrings
         if isinstance(mapping, str):
             self._fmap = fermionic_mapping[mapping]
         else:
             self._fmap = mapping
         if isinstance(self._H, MixedFermionBosonOperator):
             self._H.set_fermionic_encoding(self._fmap)
     def set_bosonic_mapping(self, mapping, max_bosons_per_state):
+        # self._bmap is never used for LinearCombinationOfPauliStrings
         if isinstance(mapping, str):
             self._bmap = bosonic_mapping[mapping](max_bosons_per_state)
         else:
@@ -63,9 +117,10 @@ def num_qubits(self):
             return self._H.n_qubits
         elif isinstance(self._H, MixedFermionBosonOperator):
             return self._H.num_qubits()
+        elif isinstance(self._H, LinearCombinationOfPauliStrings):
+            return self._H.num_qubits()
         else:
-            raise TypeError(" ".join(["Unable to determine the number of qubits for types other",
-                                      "than \"InteractionOperator\"."]))
+            raise TypeError("Unable to determine the number of qubits.")
     def get_all_pauli_strings(self, return_as="tuples"):
         # Returns all Pauli strings as a flat data structure, specifically a dictionary where the
         # key is the Pauli string and the value is the coefficient.
@@ -76,6 +131,8 @@ def get_all_pauli_strings(self, return_as="tuples"):
         # -- If return_as == "strings": The Pauli string is encoded as a character string, where
         #    each character is a Pauli matrix, explicitly including identity entries.  For example,
         #    assuming 6 qubits, "XIIZII".
+        # TODO: I'd prefer that the flag identify not the data structure but the concept: dense vs
+        #       sparse, rather than strings vs tuples.
         if return_as == "tuples":
             if isinstance(self._H, InteractionOperator):
                 return self._fmap(self._H).terms
@@ -84,14 +141,20 @@ def get_all_pauli_strings(self, return_as="tuples"):
                 #       encodings selected by Hamiltonian.  Clean this up.  Probably by deferring
                 #       the specification of encodings for MixedFermionBosonOperator?
                 return self._H.generate_qubit_operator().terms
+            elif isinstance(self._H, LinearCombinationOfPauliStrings):
+                return self._H.get_sparse_pauli_strings()
             else:
                 raise TypeError(
                     f"Unable to generate Pauli strings from object of type \"{type(self._H)}\".")
         elif return_as == "strings":
-            as_tuples = self.get_all_pauli_strings(return_as="tuples")
-            Nq = self.num_qubits()
-            return (tuple_to_string(pauli, coef, Nq)
-                    for pauli, coef in as_tuples.items() if pauli != ())
+            if isinstance(self._H, LinearCombinationOfPauliStrings):
+                return self._H.get_dense_pauli_strings()
+            else:
+                as_tuples = self.get_all_pauli_strings(return_as="tuples")
+                Nq = self.num_qubits()
+                # TODO: Why is the all-identity string excluded?  Isn't that a bug?
+                return {sparse_to_dense_pauli(pauli, Nq) : coef
+                    for pauli, coef in as_tuples.items() if pauli != ()}
         else:
             raise ValueError("  ".join([
                 "The value of return_as must be \"tuples\" or \"strings\".",
@@ -125,6 +188,8 @@ def energy_shift(self, dE):
             self._H = InteractionOperator(t0, t1, t2)
         elif isinstance(self._H, MixedFermionBosonOperator):
             self._H.energy_shift(dE)
+        elif isinstance(self._H, LinearCombinationOfPauliStrings):
+            self._H.energy_shift(dE)
         else:
             raise TypeError(
                     f"Unable to shift a fermionic Hamiltonian of type \"{type(self._H)}\".")
@@ -154,7 +219,6 @@ def compute_initial_energy_bounds(
             Ehi1 = config_hamiltonian.upper_bound
         config_general.log_verbose(f"-- initial bounds = [{Elo1}, {Ehi1})")
         return (Elo1, Ehi1)
-    # TODO: Should there be a method to generate a pyLIQTR problem instance?
 
 # -------------------------------------------------------------------------------------------------
 
@@ -236,6 +300,64 @@ def get_optional_scalar(name, default_value):
 
 # -------------------------------------------------------------------------------------------------
 
+def load_pauli(
+        config_general: GeneralConfiguration,
+        config_hamiltonian: HamiltonianConfiguration):
+    filename = config_hamiltonian.filename
+    config_general.log(
+            f"Loading Pauli string Hamiltonian from file \"{filename}\".")
+    extension = filename[filename.rfind('.')+1:]
+    if extension in [ "txt", "dat" ]:
+        fmt = None
+        numq = 0
+        pauli_dict = dict()
+        with open(filename, 'r') as file:
+            for line in file:
+                line = line.strip()
+                if line[0] == "#":
+                    continue
+                idx = line.find(' ')
+                coef_str = line[:idx].strip()
+                pauli = line[idx+1:].strip()
+                if pauli[0] == '[':
+                    if fmt is not None and fmt != "sparse":
+                        raise ValueError("Inconsistent Pauli string file format.")
+                    fmt = "sparse"
+                    coefficient = np.complex128(coef_str[1:-1])
+                    if pauli[-1] == '+':
+                        pauli = pauli[:pauli.rfind(']')+1]
+                    pauli = pauli[1:-1]
+                    pauli_tokens = pauli.split()
+                    sparse_pauli = tuple()
+                    for token in pauli_tokens:
+                        op = token[0]
+                        idx = int(token[1:])
+                        numq = max(numq, idx+1)
+                        sparse_pauli = (*sparse_pauli, (idx, op))
+                    pauli_dict[sparse_pauli] = coefficient
+                else:
+                    if fmt is not None and fmt != "dense":
+                        raise ValueError("Inconsistent Pauli string file format.")
+                    fmt = "dense"
+                    coefficient = np.complex128(coef_str)
+                    if numq != 0 and len(pauli) != numq:
+                        raise ValueError("Inconsistent dense Pauli string length.")
+                    numq = len(pauli)
+                    pauli_dict[pauli] = coefficient
+        if fmt == "dense":
+            return Hamiltonian(LinearCombinationOfPauliStrings(num_qubits=numq, dense=pauli_dict))
+        elif fmt == "sparse":
+            return Hamiltonian(LinearCombinationOfPauliStrings(num_qubits=numq, sparse=pauli_dict))
+        else:
+            raise ValueError(f"Invalid Pauli format: \"{fmt}\".")
+    elif extension == "json":
+        raise NotImplementedError("JSON Pauli string file not yet implemented.")
+    else:
+        raise ValueError(
+            f"Invalid file extension for loading a Pauli string file: \"{extension}\".")
+
+# -------------------------------------------------------------------------------------------------
+
 # TODO: Scott and I have both spent time chasing down the types of different things for a variety
 #       of reasons.  If we can get this code to the point where it always returns the same type
 #       regardless of the options passed in, then we should annotate the return type.  If it turns
@@ -255,5 +377,7 @@ def get_physical_hamiltonian(
         return load_LCPS(config_general, config_hamiltonian)
     elif config_hamiltonian.source == "hdf5":
         return load_hdf5(config_general, config_hamiltonian)
+    elif config_hamiltonian.source == "pauli":
+        return load_pauli(config_general, config_hamiltonian)
     else:
         raise ValueError(f"Invalid Hamiltonian source \"{config_hamiltonian.source}\".")

qre/qre_types.py

@@ -52,20 +52,30 @@ def __init__(self):
         self.upper_bound = float('inf')
         self.exact_energy_lower_bound = False
         self.exact_energy_upper_bound = False
+        # TODO: The F2Q and B2Q transforms should be arguments to the second-quantization file
+        #       loading function(s)
         self.fermion_to_qubit_transform = "JW"
         self.boson_to_qubit_transform = "binary"
     def _only_once(self):
         if self.source is not None:
             print("Already set Hamiltonian source to {self.source}.")
             assert self.source is None
+    # TODO: Rename to indicate this is second-quantization tensors
     def load_numpy(self, filename):
         self._only_once()
         self.source = "numpy"
         self.filename = filename
+    # TODO: Rename to indicate this is second-quantization tensors
+    # TODO: Can we merge HDF5 and NumPy paths by inspecting the extensions and looking for HDF5,
+    #       NPY, or NPY?
     def load_hdf5(self, filename):
         self._only_once()
         self.source = "hdf5"
         self.filename = filename
+    def load_pauli_strings(self, filename):
+        self._only_once()
+        self.source = "pauli"
+        self.filename = filename
     def set_energy_lower_bound(self, value, exact=False):
         self.lower_bound = value
         self.exact_energy_lower_bound = exact

qre/qre_unitary.py

@@ -120,7 +120,7 @@ def encode_ramped_trotter(
     config_general.log(f"-- using {method} Trotter formula with {Nsteps} steps ({Nsteps0})")
 
     return build_ramped_trotterized_unitary(
-            hamiltonian.get_all_pauli_strings(return_as='strings'),
+            hamiltonian.get_all_pauli_strings(return_as='strings').items(),
             method,
             timestep,
             Nsteps)