Add function to generate modulated Gaussian control

qctrlopencontrols/driven_controls/predefined.py

@@ -38,6 +38,7 @@
     WAMF1,
 )
 from ..exceptions import ArgumentsValueError
+from ..utils import check_arguments
 from .driven_control import DrivenControl
 
 
@@ -812,3 +813,128 @@ def _new_walsh_amplitude_modulated_filter_1_control(
         durations=durations,
         **kwargs,
     )
+
+
+def new_modulated_gaussian_control(
+    maximum_rabi_rate: float,
+    minimum_segment_duration: float,
+    duration: float,
+    modulation_frequency: float,
+) -> DrivenControl:
+    """
+    Generate a Gaussian driven control sequence modulated by a sinusoidal signal at a specific
+    frequency.
+
+    The net effect of this control sequence is an identity gate.
+
+    Parameters
+    ----------
+    maximum_rabi_rate: float
+        Maximum Rabi rate of the system.
+
+    minimum_segment_duration : float
+        Minimum length of each segment in the control sequence.
+
+    duration : float
+        Total duration of the control sequence.
+
+    modulation_frequency: float
+        Frequency of the modulation sinusoidal signal.
+
+    Returns
+    -------
+    DrivenControl
+        A control sequence as an instance of DrivenControl.
+    """
+
+    check_arguments(
+        maximum_rabi_rate > 0.0,
+        "Maximum Rabi rate must be greater than zero.",
+        {"maximum_rabi_rate": maximum_rabi_rate},
+    )
+
+    check_arguments(
+        minimum_segment_duration > 0.0,
+        "Minimum segment duration must be greater than zero.",
+        {"minimum_segment_duration": minimum_segment_duration},
+    )
+
+    check_arguments(
+        duration > minimum_segment_duration,
+        "Total duration must be greater than minimum segment duration.",
+        {"duration": duration, "minimum_segment_duration": minimum_segment_duration,},
+    )
+
+    # default spread of the gaussian shaped pulse as a fraction of its duration
+    _pulse_width = 0.1
+
+    # default mean of the gaussian shaped pulse as a fraction of its duration
+    _pulse_mean = 0.5
+
+    min_required_upper_bound = np.sqrt(2 * np.pi) / (_pulse_width * duration)
+    check_arguments(
+        maximum_rabi_rate >= min_required_upper_bound,
+        "Maximum Rabi rate must be large enough to permit a 2Pi rotation.",
+        {"maximum_rabi_rate": maximum_rabi_rate},
+        extras={
+            "minimum required value for upper_bound "
+            "(sqrt(2pi)/(0.1*maximum_duration))": min_required_upper_bound
+        },
+    )
+
+    # work out exact segment duration
+    segment_count = int(np.ceil(duration / minimum_segment_duration))
+    segment_duration = duration / segment_count
+    segment_start_times = np.arange(segment_count) * segment_duration
+    segment_midpoints = segment_start_times + segment_duration / 2
+
+    # prepare a base gaussian shaped pulse
+    gaussian_mean = _pulse_mean * duration
+    gaussian_width = _pulse_width * duration
+    base_gaussian_segments = np.exp(
+        -0.5 * ((segment_midpoints - gaussian_mean) / gaussian_width) ** 2
+    )
+
+    if modulation_frequency != 0:
+        # prepare the modulation signals. We use sinusoids that are zero at the center of the pulse,
+        # which ensures the pulses are antisymmetric about the center of the pulse and thus effect
+        # a net zero rotation.
+        modulation_signals = np.sin(
+            2.0 * np.pi * modulation_frequency * (segment_midpoints - duration / 2)
+        )
+        # modulate the base gaussian
+        modulated_gaussian_segments = base_gaussian_segments * modulation_signals
+
+        # maximum segment value
+        pulse_segments_maximum = np.max(modulated_gaussian_segments)
+        # normalize to maximum Rabi rate
+        modulated_gaussian_segments = (
+            maximum_rabi_rate * modulated_gaussian_segments / pulse_segments_maximum
+        )
+    else:
+        # for the zero-frequency pulse, we need to produce the largest possible full rotation (i.e.
+        # multiple of 2pi) while respecting the maximum Rabi rate. Note that if the maximum Rabi
+        # rate does not permit even a single rotation (which could happen to a small degree due to
+        # discretization issues) then we allow values to exceed the maximum Rabi rate.
+        normalized_gaussian_segments = base_gaussian_segments / np.max(
+            base_gaussian_segments
+        )
+        maximum_rotation_angle = (
+            segment_duration * np.sum(normalized_gaussian_segments) * maximum_rabi_rate
+        )
+        maximum_full_rotation_angle = max(
+            maximum_rotation_angle - maximum_rotation_angle % (2 * np.pi), 2 * np.pi
+        )
+        modulated_gaussian_segments = (
+            normalized_gaussian_segments
+            * maximum_rabi_rate
+            * (maximum_full_rotation_angle / maximum_rotation_angle)
+        )
+
+    azimuthal_angles = [0 if v >= 0 else np.pi for v in modulated_gaussian_segments]
+
+    return DrivenControl(
+        rabi_rates=np.abs(modulated_gaussian_segments),
+        azimuthal_angles=azimuthal_angles,
+        durations=np.array([segment_duration] * segment_count),
+    )

tests/test_predefined_driven_controls.py

@@ -28,9 +28,11 @@
     CORPSE_IN_SK1,
     PRIMITIVE,
     SCROFULOUS,
+    SIGMA_X,
     SK1,
     WAMF1,
 )
+from qctrlopencontrols.driven_controls.predefined import new_modulated_gaussian_control
 from qctrlopencontrols.exceptions import ArgumentsValueError
 
 
@@ -592,3 +594,117 @@ def test_walsh_control():
     )
 
     assert np.allclose(pi_on_4_segments, _pi_on_4_segments)
+
+
+def test_modulated_gaussian_control():
+    """
+    Tests modulated Gaussian control at different modulate frequencies.
+    """
+    maximum_rabi_rate = 20 * 2 * np.pi
+    minimum_segment_duration = 0.01
+    maximum_duration = 0.2
+
+    # set modulation frequency to 0
+    pulses_zero = new_modulated_gaussian_control(
+        maximum_rabi_rate=maximum_rabi_rate,
+        minimum_segment_duration=minimum_segment_duration,
+        duration=maximum_duration,
+        modulation_frequency=0,
+    )
+    pulse_zero_segments = [
+        {"duration": d, "value": np.real(v)}
+        for d, v in zip(
+            pulses_zero.durations,
+            pulses_zero.rabi_rates * np.exp(1j * pulses_zero.azimuthal_angles),
+        )
+    ]
+
+    # set modulation frequency to 50/3
+    pulses_non_zero = new_modulated_gaussian_control(
+        maximum_rabi_rate=maximum_rabi_rate,
+        minimum_segment_duration=minimum_segment_duration,
+        duration=maximum_duration,
+        modulation_frequency=50 / 3,
+    )
+    pulse_non_zero_segments = [
+        {"duration": d, "value": np.real(v)}
+        for d, v in zip(
+            pulses_non_zero.durations,
+            pulses_non_zero.rabi_rates * np.exp(1j * pulses_non_zero.azimuthal_angles),
+        )
+    ]
+
+    # pulses should have 20 segments; 0.2/0.01 = 20
+    assert len(pulse_zero_segments) == 20
+    assert len(pulse_non_zero_segments) == 20
+
+    # determine the segment mid-points
+    segment_mid_points = 0.2 / 20 * (0.5 + np.arange(20))
+    base_gaussian_mean = maximum_duration * 0.5
+    base_gaussian_width = maximum_duration * 0.1
+    base_gaussian = np.exp(
+        -0.5 * ((segment_mid_points - base_gaussian_mean) / base_gaussian_width) ** 2.0
+    ) / (np.sqrt(2 * np.pi) * base_gaussian_width)
+
+    # for modulation at frequency = 0
+    segment_values = np.array([p["value"] for p in pulse_zero_segments])
+    segment_durations = np.array([p["duration"] for p in pulse_zero_segments])
+    # The base Gaussian creates a rotation of 1rad and has maximum value
+    # 1/(sqrt(2pi)*0.2*0.1)=~19.9. Therefore, with a maximum Rabi rate of 20*2pi, we can achieve
+    # only a single 2pi rotation, which corresponds to scaling up the Gaussian by 2pi.
+    expected_gaussian = 2 * np.pi * base_gaussian
+
+    assert np.allclose(segment_values, expected_gaussian)
+    assert np.allclose(segment_durations, 0.2 / 20)
+
+    # for modulation at frequency = 50/3
+    segment_values = np.array([segment["value"] for segment in pulse_non_zero_segments])
+    segment_durations = np.array(
+        [segment["duration"] for segment in pulse_non_zero_segments]
+    )
+    expected_base_segments = base_gaussian * np.sin(
+        2 * np.pi * (50 / 3) * (segment_mid_points - 0.2 / 2)
+    )
+    expected_base_segments /= np.max(expected_base_segments)
+    expected_base_segments *= maximum_rabi_rate
+
+    assert np.allclose(segment_values, expected_base_segments)
+    assert np.allclose(segment_durations, 0.2 / 20)
+
+
+def test_modulated_gaussian_control_give_identity_gate():
+    """
+    Tests that the modulated Gaussian sequences produce identity gates.
+
+    Apply the modulated sequences to drive a noiseless qubit rotating along X. The net
+    effect should be an identity gate.
+    """
+
+    maximum_rabi_rate = 50 * 2 * np.pi
+    minimum_segment_duration = 0.02
+    maximum_duration = 0.2
+
+    pulses = [
+        new_modulated_gaussian_control(
+            maximum_rabi_rate=maximum_rabi_rate,
+            minimum_segment_duration=minimum_segment_duration,
+            duration=maximum_duration,
+            modulation_frequency=f,
+        )
+        for f in [0, 20]
+    ]
+
+    unitaries = [
+        np.linalg.multi_dot(
+            [
+                np.cos(d * v) * np.eye(2) + 1j * np.sin(d * v) * SIGMA_X
+                for d, v in zip(
+                    p.durations, p.rabi_rates * np.exp(1j * p.azimuthal_angles),
+                )
+            ]
+        )
+        for p in pulses
+    ]
+
+    for _u in unitaries:
+        assert np.allclose(_u, np.eye(2))