Support frame transformation and dressed states

src/qutip_qip/device/processor.py

@@ -36,6 +36,7 @@
 
 import numpy as np
 from scipy.interpolate import CubicSpline
+from numpy.linalg import eigh
 
 from qutip import (Qobj, QobjEvo, identity, tensor, mesolve, mcsolve)
 from ..operations import expand_operator, globalphase
@@ -50,7 +51,7 @@
 
 
 class Processor(object):
-    """
+    r"""
     A simulator of a quantum device based on the QuTiP solver
     :func:`qutip.mesolve`.
     It is defined by the available driving Hamiltonian and
@@ -97,9 +98,6 @@ class Processor(object):
     num_qubits: int
         The number of component systems.
 
-    pulses : list of :class:`.Pulse`
-        A list of control pulses of this device
-
     t1: float or list
         Characterize the decoherence of amplitude damping of
         each qubit.
@@ -115,14 +113,22 @@ class Processor(object):
     drift : :class:`.Drift`
         A `Drift` object representing the drift Hamiltonians.
 
-    dims: list
-        The dimension of each component system.
-        Default value is a
-        qubit system of ``dim=[2,2,2,...,2]``
-
-    spline_kind: str
-        Type of the coefficient interpolation.
-        See parameters of :class:`.Processor` for details.
+    use_dressed_states: bool
+        If ``True``, the drift Hamiltonian will be diagonalized
+        to obtain a dressed frame, i.e. an interaction picture, in which
+        the dressed eigenstates are used the computational basis.
+        All the Hamiltonians and collapse opeartors
+        will be transformed accordingly when generating the dynamics.
+        It will not change the saved Hamiltonian model.
+        The loaded quantum circuit and the result of evolution
+        are given with respect to the dressed states.
+        A dressed state :math:`\hat{\left|i\right\rangle}`
+        is defined as the eigenstate that
+        has the maximal overlap with the corresponding bare qubit state
+        :math:`\left| i \right\rangle`.
+        The dressed states are saved as a unitary matrix and
+        can be obtained by :obj:`.get_transformation`,
+        where each column corresponds to a dressed state.
     """
     def __init__(self, num_qubits=None, t1=None, t2=None,
                  dims=None, spline_kind="step_func", N=None):
@@ -138,6 +144,8 @@ def __init__(self, num_qubits=None, t1=None, t2=None,
             self.dims = dims
         self.pulse_mode = "discrete"
         self.spline_kind = spline_kind
+        self.use_dressed_states = False
+        self._transformation = None
 
     @property
     def N(self):
@@ -582,6 +590,55 @@ def read_coeff(self, file_name, inctime=True):
             self.coeffs = data[:, 1:].T
             return self.get_full_tlist, self.coeffs
 
+    def compute_dressed_states(self):
+        """
+        Compute the dressed state using the drift Hamiltonian.
+        The eigenstates are saved as a transformation unitary, where
+        each column corresponds to one eigenstate.
+        The transformation matrix can be obtained by
+        :obj:`.get_transformation`.
+        """
+        if self.drift is None:
+            raise ValueError("No drift Hamiltonian defined.")
+        ham = self.drift.get_ideal_qobjevo(dims=self.dims)(0)
+        U = eigh(ham.full())[1]
+        # Permute the eigenstates in U according to the overlap
+        # with bare qubit states so that the transformed Hamiltonian
+        # match the definfinition of logical qubits.
+        # A logical qubit state |i> is defined as the eigenstate that
+        # has the maximal overlap with the corresponding bare qubit state |i>.
+        qubit_state_labels = np.argmax(np.abs(U), axis=0)
+        if len(qubit_state_labels) != len(set(qubit_state_labels)):
+            raise ValueError(
+                "The definition of dressed states is ambiguous."
+                "Please define the unitary manually by the attributes"
+                "Processor.dressing_unitary")
+        U = U[:, np.argsort(qubit_state_labels)]  # Permutation
+        U = Qobj(U, dims=[self.dims, self.dims])
+        self._transformation = U
+        return self._transformation
+
+    def get_transformation(self):
+        """
+        Get the saved unitary matrix for frame transformation.
+        """
+        return self._transformation
+
+    def set_transformation(self, unitary):
+        """
+        Set a unitary transformation that will be applied to
+        the Hamiltonian and collapse operators according to the transformation
+        :math:``O_I = U^{\dagger} O U``.
+        It will not change the saved Hamiltonian model, but only apply to it
+        when generating the dynamics.
+        """
+        if self.use_dressed_states:
+            raise ValueError(
+                "To use custom transformation, "
+                "set the option use_dressed_states to False."
+                )
+        self._transformation = unitary
+
     def get_noisy_pulses(self, device_noise=False, drift=False):
         """
         It takes the pulses defined in the `Processor` and
@@ -670,6 +727,15 @@ def get_qobjevo(self, args=None, noisy=False):
             temp.append(_merge_qobjevo([c_op], full_tlist))
         c_ops = temp
 
+        # compute the dressed Hamiltonian and collapse operator
+        if self.use_dressed_states:
+            self.compute_dressed_states()
+        if self._transformation is not None:
+            U = self._transformation
+            U_dag = self._transformation.dag()
+            final_qu = U_dag * final_qu * U
+            c_ops = [U_dag * op * U for op in c_ops]
+
         if noisy:
             return final_qu, c_ops
         else:

tests/test_processor.py

@@ -32,13 +32,15 @@
 ###############################################################################
 import os
 
+import pytest
 from numpy.testing import (
     assert_, run_module_suite, assert_allclose, assert_equal)
 import numpy as np
 
 from qutip_qip.device.processor import Processor
-from qutip import (basis, sigmaz, sigmax, sigmay, identity, destroy, tensor,
-                Options, rand_ket, rand_dm, fidelity)
+from qutip import (
+    basis, sigmaz, sigmax, sigmay, identity, destroy, tensor,
+    Options, rand_ket, rand_dm, fidelity, Qobj)
 from qutip_qip.operations import hadamard_transform
 from qutip_qip.noise import (
     DecoherenceNoise, RandomNoise, ControlAmpNoise)
@@ -374,3 +376,55 @@ def test_repeated_use_of_processor(self):
         result1 = processor.run_state(basis(2, 0), tlist=np.linspace(0, 1, 10))
         result2 = processor.run_state(basis(2, 0), tlist=np.linspace(0, 1, 10))
         assert_allclose(result1.states[-1].full(), result2.states[-1].full())
+
+    def test_frame_transformation(self):
+        # Define the drift, the bare eigenvalue and the transformation U.
+        H_d = sigmaz() + 0.1 * sigmay()
+        energy, U = np.linalg.eigh(H_d.full())
+        U = Qobj(np.flip(U, axis=1))  # largest eigenvalue first
+        H_c = sigmay()
+
+        # Initialize a processor
+        processor = Processor(1)
+        processor.use_dressed_states = True
+        processor.add_drift(H_d, targets=0)
+
+        # only drift
+        processor.add_pulse(  # add a dummy empty pulse
+            Pulse(
+                0. * H_c, targets=0,
+                coeff=True, tlist=np.array([0., 1.]))
+            )
+        processor.compute_dressed_states()
+        np.testing.assert_allclose(
+            processor.get_transformation(),
+            U
+        )
+        qu, _ = processor.get_qobjevo(noisy=True)
+        np.testing.assert_allclose(
+            qu(0).full(),
+            np.diag([max(energy), min(energy)])
+        )
+
+        # drift + control
+        processor.add_pulse(
+            Pulse(
+                H_c, targets=0,
+                coeff=True, tlist=np.array([0., 1.]))
+            )
+        qu, _ = processor.get_qobjevo(noisy=True)
+        np.testing.assert_allclose(
+            qu(0).full(),
+            (U.dag() * (H_d + H_c) * U).full()
+        )
+
+        # set custom transformation
+        processor.use_dressed_states = False
+        processor.set_transformation(sigmax())
+        assert(processor.get_transformation() == sigmax())
+
+        # Ambiguous dressed states
+        processor = Processor(1)
+        processor.add_drift(sigmax(), targets=0)
+        with pytest.raises(ValueError):
+            processor.compute_dressed_states()