Migrating Marek's code for DBI to discontinue the old repository

src/qibo/models/dbi/double_bracket.py

@@ -5,7 +5,6 @@
 import hyperopt
 import numpy as np
 
-from qibo.config import raise_error
 from qibo.hamiltonians import Hamiltonian
 
 
@@ -33,13 +32,12 @@ class DoubleBracketIteration:
     Example:
         .. testcode::
 
-            import numpy as np
             from qibo.models.dbi.double_bracket import DoubleBracketIteration, DoubleBracketGeneratorType
-            from qibo.hamiltonians import Hamiltonian
             from qibo.quantum_info import random_hermitian
+            from qibo.hamiltonians import Hamiltonian
 
             nqubits = 4
-            h0 = random_hermitian(2**nqubits)
+            h0 = random_hermitian(2**nqubits, seed=2)
             dbf = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0))
 
             # diagonalized matrix
@@ -68,14 +66,14 @@ def __call__(
             )
         elif mode is DoubleBracketGeneratorType.single_commutator:
             if d is None:
-                raise_error(ValueError, f"Cannot use group_commutator with matrix {d}")
+                d = self.diagonal_h_matrix
             operator = self.backend.calculate_matrix_exp(
                 1.0j * step,
                 self.commutator(d, self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.group_commutator:
             if d is None:
-                raise_error(ValueError, f"Cannot use group_commutator with matrix {d}")
+                d = self.diagonal_h_matrix
             operator = (
                 self.h.exp(-step)
                 @ self.backend.calculate_matrix_exp(-step, d)
@@ -103,12 +101,12 @@ def off_diag_h(self):
 
     @property
     def off_diagonal_norm(self):
-        """Norm of off-diagonal part of H matrix."""
+        r"""Hilbert Schmidt norm of off-diagonal part of H matrix, namely :math:`\\text{Tr}(\\sqrt{A^{\\dagger} A})`."""
         off_diag_h_dag = self.backend.cast(
             np.matrix(self.backend.to_numpy(self.off_diag_h)).getH()
         )
-        return np.real(
-            np.trace(self.backend.to_numpy(off_diag_h_dag @ self.off_diag_h))
+        return np.sqrt(
+            np.real(np.trace(self.backend.to_numpy(off_diag_h_dag @ self.off_diag_h)))
         )
 
     @property
@@ -125,6 +123,7 @@ def hyperopt_step(
         optimizer: callable = None,
         look_ahead: int = 1,
         verbose: bool = False,
+        d: np.array = None,
     ):
         """
         Optimize iteration step.
@@ -136,7 +135,8 @@ def hyperopt_step(
             space: see hyperopt.hp possibilities;
             optimizer: see hyperopt algorithms;
             look_ahead: number of iteration steps to compute the loss function;
-            verbose: level of verbosity.
+            verbose: level of verbosity;
+            d: diagonal operator for generating double-bracket iterations.
 
         Returns:
             (float): optimized best iteration step.
@@ -148,28 +148,28 @@ def hyperopt_step(
 
         space = space("step", step_min, step_max)
         best = hyperopt.fmin(
-            fn=partial(self.loss, look_ahead=look_ahead),
+            fn=partial(self.loss, d=d, look_ahead=look_ahead),
             space=space,
             algo=optimizer.suggest,
             max_evals=max_evals,
             verbose=verbose,
         )
-
         return best["step"]
 
-    def loss(self, step: float, look_ahead: int = 1):
+    def loss(self, step: float, d: np.array = None, look_ahead: int = 1):
         """
         Compute loss function distance between `look_ahead` steps.
 
         Args:
             step: iteration step.
+            d: diagonal operator, use canonical by default.
             look_ahead: number of iteration steps to compute the loss function;
         """
         # copy initial hamiltonian
         h_copy = deepcopy(self.h)
 
         for _ in range(look_ahead):
-            self.__call__(mode=self.mode, step=step)
+            self.__call__(mode=self.mode, step=step, d=d)
 
         # off_diagonal_norm's value after the steps
         loss = self.off_diagonal_norm

src/qibo/models/dbi/utils.py

@@ -0,0 +1,159 @@
+from copy import deepcopy
+from itertools import product
+from typing import Optional
+
+import numpy as np
+from hyperopt import hp, tpe
+
+from qibo import symbols
+from qibo.hamiltonians import SymbolicHamiltonian
+from qibo.models.dbi.double_bracket import (
+    DoubleBracketGeneratorType,
+    DoubleBracketIteration,
+)
+
+
+def generate_Z_operators(nqubits: int):
+    """Generate a dictionary containing 1) all possible products of Pauli Z operators for L = n_qubits and 2) their respective names.
+    Return: Dictionary with operator names (str) as keys and operators (np.array) as values
+
+     Example:
+        .. testcode::
+
+            from qibo.models.dbi.utils import generate_Z_operators
+            from qibo.models.dbi.double_bracket import DoubleBracketIteration
+            from qibo.quantum_info import random_hermitian
+            from qibo.hamiltonians import Hamiltonian
+            import numpy as np
+
+            nqubits = 4
+            h0 = random_hermitian(2**nqubits)
+            dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0))
+            generate_Z = generate_Z_operators(nqubits)
+            Z_ops = list(generate_Z.values())
+
+            delta_h0 = dbi.diagonal_h_matrix
+            dephasing_channel = (sum([Z_op @ h0 @ Z_op for Z_op in Z_ops])+h0)/2**nqubits
+            norm_diff = np.linalg.norm(delta_h0 - dephasing_channel)
+    """
+    # list of tuples, e.g. ('Z','I','Z')
+    combination_strings = product("ZI", repeat=nqubits)
+    output_dict = {}
+
+    for zi_string_combination in combination_strings:
+        # except for the identity
+        if "Z" in zi_string_combination:
+            op_name = "".join(zi_string_combination)
+            tensor_op = str_to_symbolic(op_name)
+            # append in output_dict
+            output_dict[op_name] = SymbolicHamiltonian(tensor_op).dense.matrix
+    return output_dict
+
+
+def str_to_symbolic(name: str):
+    """Convert string into symbolic hamiltonian.
+    Example:
+        .. testcode::
+
+            from qibo.models.dbi.utils import str_to_symbolic
+            op_name = "ZYXZI"
+            # returns 5-qubit symbolic hamiltonian
+            ZIXZI_op = str_to_symbolic(op_name)
+    """
+    tensor_op = 1
+    for qubit, char in enumerate(name):
+        tensor_op *= getattr(symbols, char)(qubit)
+    return tensor_op
+
+
+def select_best_dbr_generator(
+    dbi_object: DoubleBracketIteration,
+    d_list: list,
+    step: Optional[float] = None,
+    step_min: float = 1e-5,
+    step_max: float = 1,
+    max_evals: int = 200,
+    compare_canonical: bool = True,
+):
+    """Selects the best double bracket rotation generator from a list and execute the rotation.
+
+    Args:
+        dbi_object (`DoubleBracketIteration`): the target DoubleBracketIteration object.
+        d_list (list): list of diagonal operators (np.array) to run from.
+        step (float): fixed iteration duration.
+            Defaults to ``None``, uses hyperopt.
+        step_min (float): minimally allowed iteration duration.
+        step_max (float): maximally allowed iteration duration.
+        max_evals (int): maximally allowed number of evaluation in hyperopt.
+        compare_canonical (bool): if `True`, the optimal diagonal operator chosen from "d_list" is compared with the canonical bracket.
+
+    Returns:
+        The updated dbi_object, index of the optimal diagonal operator, respective step duration, and evolution direction.
+    """
+    norms_off_diagonal_restriction = [
+        dbi_object.off_diagonal_norm for _ in range(len(d_list))
+    ]
+    optimal_steps, flip_list = [], []
+    for i, d in enumerate(d_list):
+        # prescribed step durations
+        dbi_eval = deepcopy(dbi_object)
+        flip_list.append(cs_angle_sgn(dbi_eval, d))
+        if flip_list[i] != 0:
+            if step is None:
+                step_best = dbi_eval.hyperopt_step(
+                    d=flip_list[i] * d,
+                    step_min=step_min,
+                    step_max=step_max,
+                    space=hp.uniform,
+                    optimizer=tpe,
+                    max_evals=max_evals,
+                )
+            else:
+                step_best = step
+            dbi_eval(step=step_best, d=flip_list[i] * d)
+            optimal_steps.append(step_best)
+            norms_off_diagonal_restriction[i] = dbi_eval.off_diagonal_norm
+    # canonical
+    if compare_canonical is True:
+        flip_list.append(1)
+        dbi_eval = deepcopy(dbi_object)
+        dbi_eval.mode = DoubleBracketGeneratorType.canonical
+        if step is None:
+            step_best = dbi_eval.hyperopt_step(
+                step_min=step_min,
+                step_max=step_max,
+                space=hp.uniform,
+                optimizer=tpe,
+                max_evals=max_evals,
+            )
+        else:
+            step_best = step
+        dbi_eval(step=step_best)
+        optimal_steps.append(step_best)
+        norms_off_diagonal_restriction.append(dbi_eval.off_diagonal_norm)
+    # find best d
+    idx_max_loss = np.argmin(norms_off_diagonal_restriction)
+    flip = flip_list[idx_max_loss]
+    step_optimal = optimal_steps[idx_max_loss]
+    dbi_eval = deepcopy(dbi_object)
+    if idx_max_loss == len(d_list) and compare_canonical is True:
+        # canonical
+        dbi_eval(step=step_optimal, mode=DoubleBracketGeneratorType.canonical)
+
+    else:
+        d_optimal = flip * d_list[idx_max_loss]
+        dbi_eval(step=step_optimal, d=d_optimal)
+    return dbi_eval, idx_max_loss, step_optimal, flip
+
+
+def cs_angle_sgn(dbi_object, d):
+    """Calculates the sign of Cauchy-Schwarz Angle :math:`\\langle W(Z), W({\\rm canonical}) \\rangle_{\\rm HS}`."""
+    norm = np.trace(
+        np.dot(
+            np.conjugate(
+                dbi_object.commutator(dbi_object.diagonal_h_matrix, dbi_object.h.matrix)
+            ).T,
+            dbi_object.commutator(d, dbi_object.h.matrix),
+        )
+    )
+    return np.sign(norm)