Add support for time step dt != 1 ns

sequencing/modes.py

@@ -607,6 +607,17 @@ def initialize(self):
         super().initialize()
         self.add_pulse(cls=SmoothedConstantPulse)
         self.add_pulse(cls=GaussianPulse)
+        self._dt = 1
+
+    @property
+    def dt(self):
+        return self._dt
+
+    @dt.setter
+    def dt(self, dt):
+        self._dt = dt
+        for pulse in self.pulses.values():
+            pulse.dt = dt
 
     def add_pulse(self, cls=GaussianPulse, name=None, error_if_exists=False, **kwargs):
         """Creates a new pulse of type ``cls`` and adds it to ``self.pulses``.
@@ -739,7 +750,7 @@ def rotate(self, angle, phase, pulse_name=None, unitary=False, **kwargs):
             f"{self.name}.x": HTerm(self.x, c_wave.real),
             f"{self.name}.y": HTerm(self.y, c_wave.imag),
         }
-        return Operation(len(c_wave), terms)
+        return Operation(len(c_wave) * self.dt, terms)
 
     def rotate_x(self, angle, unitary=False, **kwargs):
         """Generate a rotation about the x axis.
@@ -883,4 +894,4 @@ def displace(self, alpha, unitary=False, pulse_name=None, **kwargs):
             f"{self.name}.x": HTerm(self.x, i),
             f"{self.name}.y": HTerm(self.y, q),
         }
-        return Operation(len(c_wave), terms)
+        return Operation(len(c_wave) * self.dt, terms)

sequencing/pulses.py

@@ -34,6 +34,7 @@ def array_pulse(
     amp=1,
     phase=0,
     detune=0,
+    dt=1.0,
     noise_sigma=0,
     noise_alpha=0,
     scale_noise=False,
@@ -47,9 +48,10 @@ def array_pulse(
             If None, the imaginary part is set to 0. Default: None.
         amp (float): Factor by which to scale the waveform amplitude.
             Default: 1.
-        phase (optionla, float): Phase offset in radians. Default: 0.
+        phase (optional, float): Phase offset in radians. Default: 0.
         detune (optional, float): Software detuning/time-dependent phase to
             apply to the waveform, in GHz. Default: 0.
+        dt (optional, float): Sample time in nanoseconds.
         noise_sigma (optional, float): Standard deviation of additive Gaussian
             noise applied to the pulse (see scale_noise).
             Default: 0.
@@ -68,7 +70,7 @@ def array_pulse(
     if q_wave is None:
         q_wave = np.zeros_like(i_wave)
     if detune:
-        ts = np.arange(len(i_wave))
+        ts = np.linspace(0, len(i_wave) * dt, len(i_wave))
         c_wave = (i_wave + 1j * q_wave) * np.exp(-2j * np.pi * ts * detune)
         i_wave, q_wave = c_wave.real, c_wave.imag
     if noise_sigma:
@@ -87,15 +89,17 @@ def array_pulse(
     return c_wave
 
 
-def gaussian_wave(sigma, chop=4):
-    ts = np.linspace(-chop // 2 * sigma, chop // 2 * sigma, int(chop * sigma // 4) * 4)
+def gaussian_wave(sigma, chop=4, dt=1):
+    length = chop * sigma
+    ts = np.linspace(-length / 2, length / 2, int(length / dt))
     P = np.exp(-(ts**2) / (2.0 * sigma**2))
     ofs = P[0]
     return (P - ofs) / (1 - ofs)
 
 
-def gaussian_deriv_wave(sigma, chop=4):
-    ts = np.linspace(-chop // 2 * sigma, chop // 2 * sigma, int(chop * sigma // 4) * 4)
+def gaussian_deriv_wave(sigma, chop=4, dt=1):
+    length = chop * sigma
+    ts = np.linspace(-length / 2, length / 2, int(length / dt))
     ofs = np.exp(-ts[0] ** 2 / (2 * sigma**2))
     return (0.25 / sigma**2) * ts * np.exp(-(ts**2) / (2 * sigma**2)) / (1 - ofs)
 
@@ -128,96 +132,98 @@ def ring_up_wave(length, reverse=False, shape="tanh", **kwargs):
     return wave
 
 
-def ring_up_gaussian_flattop(length, sigma, ramp_offset=None):
+def ring_up_gaussian_flattop(length, sigma, ramp_offset=None, dt=1):
     ramp_offset = 0 if ramp_offset is None else ramp_offset
 
-    def _ring_up(ts):
-        if np.abs(ts) < ramp_offset:
+    def _ring_up(t):
+        if np.abs(t) < ramp_offset:
             return 1.0
         elif ts > ramp_offset:
-            return np.exp(-((ts - ramp_offset) ** 2) / (2.0 * sigma**2))
-        else:  # ts < ramp_offset
-            return np.exp(-((ts + ramp_offset) ** 2) / (2.0 * sigma**2))
+            return np.exp(-((t - ramp_offset) ** 2) / (2.0 * sigma**2))
+        else:  # t < ramp_offset
+            return np.exp(-((t + ramp_offset) ** 2) / (2.0 * sigma**2))
 
-    ts = np.linspace(-length + 1, 0, length)
+    ts = np.linspace(-length + 1, 0, int(length / dt))
     P = np.array([_ring_up(t) for t in ts])
     # normalize so tail amp = 0 and max amp = 0
     ofs = P[0]
     return (P - ofs) / (1 - ofs)
 
 
-def ring_up_cos(length):
-    return 0.5 * (1 - np.cos(np.linspace(0, np.pi, length)))
+def ring_up_cos(length, dt=1):
+    return 0.5 * (1 - np.cos(np.linspace(0, np.pi, int(length / dt))))
 
 
-def ring_up_tanh(length):
-    ts = np.linspace(-2, 2, length)
+def ring_up_tanh(length, dt=1):
+    ts = np.linspace(-2, 2, int(length / dt))
     return (1 + np.tanh(ts)) / 2
 
 
 def smoothed_constant_wave(length, sigma, shape="tanh", **kwargs):
+    dt = kwargs.get("dt", 1)
     if sigma == 0:
-        return np.ones(length)
+        return np.ones(int(length / dt))
 
     return np.concatenate(
         [
             ring_up_wave(sigma, shape=shape, **kwargs),
-            np.ones(length - 2 * sigma),
+            np.ones(int((length - 2 * sigma) / dt)),
             ring_up_wave(sigma, reverse=True, shape=shape, **kwargs),
         ]
     )
 
 
-def constant_pulse(length=None):
+def constant_pulse(length=None, dt=1):
+    length = int(length / dt)
     i_wave, q_wave = np.ones(length), np.zeros(length)
     return i_wave, q_wave
 
 
-def gaussian_pulse(sigma=None, chop=4, drag=0):
-    i_wave = gaussian_wave(sigma, chop=chop)
-    q_wave = drag * gaussian_deriv_wave(sigma, chop=chop)
+def gaussian_pulse(sigma=None, chop=4, drag=0, dt=1):
+    i_wave = gaussian_wave(sigma, chop=chop, dt=dt)
+    q_wave = drag * gaussian_deriv_wave(sigma, chop=chop, dt=dt)
     return i_wave, q_wave
 
 
-def smoothed_constant_pulse(length=None, sigma=None, shape="tanh"):
-    i_wave = smoothed_constant_wave(length, sigma, shape=shape)
+def smoothed_constant_pulse(length=None, sigma=None, shape="tanh", dt=1):
+    i_wave = smoothed_constant_wave(length, sigma, shape=shape, dt=dt)
     q_wave = np.zeros_like(i_wave)
     return i_wave, q_wave
 
 
-def sech_wave(sigma, chop=4):
+def sech_wave(sigma, chop=4, dt=1):
     # https://arxiv.org/pdf/1704.00803.pdf
     # https://doi.org/10.1103/PhysRevA.96.042339
     rho = np.pi / (2 * sigma)
     t0 = chop * sigma // 2
-    ts = np.linspace(-t0, t0, int(chop * sigma // 4) * 4)
+    ts = np.linspace(-t0, t0, int(chop * sigma / dt))
     P = 1 / np.cosh(rho * ts)
     ofs = P[0]
     return (P - ofs) / (1 - ofs)
 
 
-def sech_deriv_wave(sigma, chop=4):
+def sech_deriv_wave(sigma, chop=4, dt=1):
     rho = np.pi / (2 * sigma)
     t0 = chop * sigma // 2
-    ts = np.linspace(-t0, t0, int(chop * sigma // 4) * 4)
+    ts = np.linspace(-t0, t0, int(chop * sigma / dt))
     ofs = 1 / np.cosh(rho * ts[0])
     P = -np.sinh(rho * ts) / np.cosh(rho * ts) ** 2
     return (P - ofs) / (1 - ofs)
 
 
-def sech_pulse(sigma=None, chop=4, drag=0):
-    i_wave = sech_wave(sigma, chop=chop)
+def sech_pulse(sigma=None, chop=4, drag=0, dt=1):
+    i_wave = sech_wave(sigma, chop=chop, dt=dt)
     # q_wave = drag * sech_deriv_wave(sigma, chop=chop)
-    q_wave = drag * np.gradient(i_wave)
+    q_wave = drag * np.gradient(i_wave) / dt
     return i_wave, q_wave
 
 
-def slepian_pulse(tau=None, width=10, drag=0):
+def slepian_pulse(tau=None, width=10, drag=0, dt=1):
     # bandwidth is relative, i.e. scaled by 1/tau
     from scipy.signal.windows import slepian
 
-    i_wave = slepian(tau, width / tau)
-    q_wave = drag * np.gradient(i_wave)
+    i_wave = slepian(tau, width / tau / dt)
+    q_wave = drag * np.gradient(i_wave) / dt
     return i_wave, q_wave
 
 
@@ -248,6 +254,7 @@ class Pulse(Parameterized):
     amp = FloatParameter(1)
     detune = GigahertzParameter(0)
     phase = RadianParameter(0)
+    dt = NanosecondParameter(1)
     noise_sigma = FloatParameter(0)
     noise_alpha = FloatParameter(0)
     scale_noise = BoolParameter(False)
@@ -295,8 +302,10 @@ def plot(self, ax=None, grid=True, legend=True, **kwargs):
         if ax is None:
             _, ax = plt.subplots()
         c_wave = self(**kwargs)
-        (line,) = ax.plot(c_wave.real, ls="-", label=self.name)
-        ax.plot(c_wave.imag, color=line._color, ls="--")
+        dt = kwargs.get("dt", self.dt)
+        ts = np.linspace(0, len(c_wave) * dt, len(c_wave))
+        (line,) = ax.plot(ts, c_wave.real, ls="-", label=self.name)
+        ax.plot(ts, c_wave.imag, color=line._color, ls="--")
         ax.grid(grid)
         if legend:
             ax.legend(loc="best")

sequencing/sequencing/basic.py

@@ -24,7 +24,7 @@ class HamiltonianChannels(object):
             collapse operators instead of Hamiltonian terms. Default: None.
     """
 
-    def __init__(self, channels=None, collapse_channels=None, t0=0):
+    def __init__(self, channels=None, collapse_channels=None, t0=0, dt=1):
         self.channels = {}
         if channels is None:
             channels = {}
@@ -41,7 +41,7 @@ def __init__(self, channels=None, collapse_channels=None, t0=0):
             self.add_channel(name, C_op=op, time_dependent=time_dependent)
         self.tmin = t0
         self.tmax = t0
-        self.dt = 1
+        self.dt = dt
 
     @property
     def times(self):
@@ -322,9 +322,13 @@ def __init__(self, system=None, channels=None, t0=0):
     def set_system(self, system, channels=None, t0=0):
         self.system = system
         self.modes = system.modes
-        self.hc = HamiltonianChannels(channels=channels, t0=t0)
+        self.hc = HamiltonianChannels(channels=channels, t0=t0, dt=system.dt)
         self._t = self.hc.tmax
 
+    @property
+    def dt(self):
+        return self.hc.dt
+
     @property
     def times(self):
         return self.hc.times

sequencing/system.py

@@ -120,9 +120,20 @@ def initialize(self):
         super().initialize()
         if self.order_modes:
             self.modes = sort_modes(self.modes)
+        self._dt = 1
         self.active_modes = self.modes
         self.coupling_terms = defaultdict(list)
 
+    @property
+    def dt(self):
+        return self._dt
+
+    @dt.setter
+    def dt(self, dt):
+        self._dt = dt
+        for mode in self.modes:
+            mode.dt = dt
+
     def __getattribute__(self, name):
         # Access modes like system.qubit.
         # System.modes can be changed at any time, so we cannot

sequencing/test/test_calibration.py

@@ -16,6 +16,7 @@ def tearDownClass(cls):
     def test_rabi_two_levels(self):
         qubit = Transmon("qubit", levels=2)
         system = System("system", modes=[qubit])
+        system.dt = 0.75
         for _ in range(5):
             _, old_amp, new_amp = tune_rabi(system, qubit.fock(0))
         self.assertLess(abs(old_amp - new_amp), 1e-7)

sequencing/test/test_modes.py

@@ -203,6 +203,10 @@ def test_rotations_pulse(self):
         q1.gaussian_pulse.sigma = 40
         system = System("system", modes=[q0, q1])
         init_state = system.fock()
+        dt = 0.5
+        system.dt = dt
+        self.assertEqual(q0.dt, dt)
+        self.assertEqual(q1.dt, dt)
 
         for qubit in [q0, q1]:
             for _ in range(1):
@@ -367,6 +371,10 @@ def test_displacement(self):
         c0 = Cavity("c0", levels=10, kerr=-10e-6)
         c1 = Cavity("c1", levels=12, kerr=-10e-6)
         system = System("system", modes=[c0, c1])
+        dt = 2
+        system.dt = dt
+        self.assertEqual(c0.dt, dt)
+        self.assertEqual(c1.dt, dt)
         init_state = system.fock()
         for cavity in [c0, c1]:
             for _ in range(1):