Add option to use QobjEvo and

sequencing/calibration.py

@@ -122,7 +122,7 @@ def tune_rabi(
         seq = get_sequence(system)
         with qubit.pulse_scale(amp, pulse_name=pulse_name):
             qubit.rotate_x(np.pi, pulse_name=pulse_name)
-        result = seq.run(init_state, e_ops=e_ops)
+        result = seq.run(init_state, e_ops=e_ops, only_final_state=False)
         e_pop.append(result.expect[0][-1])
 
     fit_result = fit_sine(amps, e_pop)
@@ -232,14 +232,14 @@ def tune_drag(
         qubit.rotate_x(np.pi, pulse_name=pulse_name)
         sync()
         qubit.rotate_y(np.pi / 2, pulse_name=pulse_name)
-        result = seq.run(init_state, e_ops=e_ops)
+        result = seq.run(init_state, e_ops=e_ops, only_final_state=False)
         XpY9.append(result.expect[0][-1])
 
         seq = get_sequence(system)
         qubit.rotate_y(np.pi, pulse_name=pulse_name)
         sync()
         qubit.rotate_x(np.pi / 2, pulse_name=pulse_name)
-        result = seq.run(init_state, e_ops=e_ops)
+        result = seq.run(init_state, e_ops=e_ops, only_final_state=False)
         YpX9.append(result.expect[0][-1])
 
     r0 = fit_line(drags, XpY9)
@@ -338,7 +338,7 @@ def tune_displacement(
         seq = get_sequence(system)
         with cavity.pulse_scale(amp):
             cavity.displace(1)
-        result = seq.run(init_state, e_ops=e_ops)
+        result = seq.run(init_state, e_ops=e_ops, only_final_state=False)
         zero_pop.append(result.expect[0][-1])
 
     fit_result = fit_displacement(amps, zero_pop)

sequencing/sequencing/basic.py

@@ -444,6 +444,7 @@ def run(
         c_ops=None,
         e_ops=None,
         options=None,
+        only_final_state=False,
         progress_bar=None,
     ):
         """Simulate the sequence using qutip.mesolve.
@@ -456,9 +457,13 @@ def run(
                 Default: None.
             options (optional, qutip.Options): qutip solver options.
                 Note: defaults to max_step = 1.
+            only_final_state (optional, bool): Whether to query the system's state
+                at only the initial and final times rather than at every time step
+                in the Hamiltonian. Default: False.
             progress_bar (optional, None): Whether to use qutip's progress bar.
                 Default: None (no progress bar).
 
+
         Returns:
             ``qutip.solver.Result``: qutip.Solver.Result instance.
         """
@@ -486,10 +491,16 @@ def run(
         c_ops.extend(C_ops)
         e_ops = ops2dms(e_ops)
 
+        if only_final_state:
+            mesolve_times = [times.min(), times.max()]
+            H = qutip.QobjEvo(H, tlist=times, state0=init_state, e_ops=e_ops)
+        else:
+            mesolve_times = times
+
         return qutip.mesolve(
             H,
             init_state,
-            times,
+            mesolve_times,
             c_ops,
             e_ops,
             options=options,
@@ -541,6 +552,7 @@ def propagator(
         if len(H) == 1 and not isinstance(H[0], list):
             # This case breaks qutip.propagator().
             H = H[0]
+
         props = qutip.propagator(
             H,
             times,

sequencing/sequencing/main.py

@@ -104,6 +104,7 @@ def run(
         c_ops=None,
         e_ops=None,
         options=None,
+        only_final_state=True,
         progress_bar=None,
     ):
         """Simulate the sequence using qutip.mesolve.
@@ -116,6 +117,9 @@ def run(
                 Default: None.
             options (optional, qutip.Options): qutip solver options.
                 Note: defaults to max_step = 1.
+            only_final_state (optional, bool): Whether to query the system's state
+                at only the initial and final times rather than at every time step
+                in the Hamiltonian. Default: True.
             progress_bar (optional, None): Whether to use qutip's progress bar.
                 Default: None (no progress bar).
 
@@ -127,6 +131,7 @@ def run(
             c_ops=c_ops,
             e_ops=e_ops,
             options=options,
+            only_final_state=only_final_state,
             progress_bar=progress_bar,
         )
 
@@ -283,7 +288,9 @@ def run(
                 item.t0 = self._t
                 seq = item.compile()
                 seq.sync()
-                result = seq.run(states[-1], options=options)
+                result = seq.run(
+                    states[-1], options=options, only_final_state=(not full_evolution)
+                )
                 if full_evolution:
                     new_states = result.states
                     new_times = result.times
@@ -295,7 +302,9 @@ def run(
                 seq = CompiledPulseSequence(self.system, t0=self._t)
                 seq.add_operation(item)
                 seq.sync()
-                result = seq.run(states[-1], options=options)
+                result = seq.run(
+                    states[-1], options=options, only_final_state=(not full_evolution)
+                )
                 if full_evolution:
                     new_states = result.states
                     new_times = result.times

sequencing/test/test_rotations.py

@@ -1,28 +1,7 @@
 import unittest
 import numpy as np
 import qutip
-from sequencing import (
-    Transmon,
-    Cavity,
-    System,
-    Sequence,
-    CTerm,
-    Operation,
-    capture_operation,
-    get_sequence,
-    sync,
-    delay,
-    delay_channels,
-    ket2dm,
-    ops2dms,
-)
-from sequencing.sequencing import (
-    ValidatedList,
-    CompiledPulseSequence,
-    PulseSequence,
-    SyncOperation,
-    HamiltonianChannels,
-)
+from sequencing import Transmon, System, Sequence, sync
 from sequencing.calibration import tune_rabi
 
 

sequencing/test/test_sequencing.py

@@ -234,7 +234,9 @@ def fit_exp_decay(xs, ys):
             return amp, tau
 
         normal_t1_result = t1_sequence(system, qubit).run(
-            qubit.fock(0), e_ops=[qubit.fock(1)]
+            qubit.fock(0),
+            e_ops=[qubit.fock(1)],
+            only_final_state=False,
         )
         t0 = qubit.gaussian_pulse.sigma * qubit.gaussian_pulse.chop
         ts = np.arange(len(normal_t1_result.states))
@@ -245,7 +247,9 @@ def fit_exp_decay(xs, ys):
         for factor in [2, 3, 4, 5]:
             pulsed_t1 = qubit.t1 / factor
             result = t1_sequence(system, qubit, pulsed_t1=pulsed_t1).run(
-                qubit.fock(0), e_ops=[qubit.fock(1)]
+                qubit.fock(0),
+                e_ops=[qubit.fock(1)],
+                only_final_state=False,
             )
             ts = np.arange(len(result.states))
             _, fit_t1 = fit_exp_decay(ts[t0:], result.expect[0][t0:])