q-optimize/c3

View on GitHub
docs/two_qubits.rst

Summary

Maintainability
Test Coverage
.. _setup-example:

Setup of a two-qubit chip with :math:`C^3`
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In this example we will set-up a two qubit quantum processor and define
a simple gate.

Imports
^^^^^^^

.. code-block:: python

    # System imports
    import copy
    import numpy as np
    import time
    import itertools
    import matplotlib.pyplot as plt
    import tensorflow as tf
    import tensorflow_probability as tfp

    # Main C3 objects
    from c3.c3objs import Quantity as Qty
    from c3.parametermap import ParameterMap as PMap
    from c3.experiment import Experiment as Exp
    from c3.model import Model as Mdl
    from c3.generator.generator import Generator as Gnr

    # Building blocks
    import c3.generator.devices as devices
    import c3.signal.gates as gates
    import c3.libraries.chip as chip
    import c3.signal.pulse as pulse
    import c3.libraries.tasks as tasks

    # Libs and helpers
    import c3.libraries.algorithms as algorithms
    import c3.libraries.hamiltonians as hamiltonians
    import c3.libraries.fidelities as fidelities
    import c3.libraries.envelopes as envelopes
    import c3.utils.qt_utils as qt_utils
    import c3.utils.tf_utils as tf_utils

Model components
^^^^^^^^^^^^^^^^

We first create a qubit. Each parameter is a Quantity (``Qty()``) object
with bounds and a unit. In :math:`C^3`, the default multi-level qubit is
a Transmon modelled as a Duffing oscillator with frequency
:math:`\omega` and anharmonicity :math:`\delta` :

.. math::

    H/\hbar = \omega b^\dagger b - \frac{\delta}{2}                        \left(b^\dagger b - 1\right) b^\dagger b

The “name” will be used to identify this qubit (or other component)
later and should thus be chosen carefully.

.. code-block:: python

    qubit_lvls = 3
    freq_q1 = 5e9
    anhar_q1 = -210e6
    t1_q1 = 27e-6
    t2star_q1 = 39e-6
    qubit_temp = 50e-3

    q1 = chip.Qubit(
        name="Q1",
        desc="Qubit 1",
        freq=Qty(
            value=freq_q1,
            min_val=4.995e9 ,
            max_val=5.005e9 ,
            unit='Hz 2pi'
        ),
        anhar=Qty(
            value=anhar_q1,
            min_val=-380e6 ,
            max_val=-120e6 ,
            unit='Hz 2pi'
        ),
        hilbert_dim=qubit_lvls,
        t1=Qty(
            value=t1_q1,
            min_val=1e-6,
            max_val=90e-6,
            unit='s'
        ),
        t2star=Qty(
            value=t2star_q1,
            min_val=10e-6,
            max_val=90e-3,
            unit='s'
        ),
        temp=Qty(
            value=qubit_temp,
            min_val=0.0,
            max_val=0.12,
            unit='K'
        )
    )

And the same for a second qubit.

.. code-block:: python

    freq_q2 = 5.6e9
    anhar_q2 = -240e6
    t1_q2 = 23e-6
    t2star_q2 = 31e-6
    q2 = chip.Qubit(
        name="Q2",
        desc="Qubit 2",
        freq=Qty(
            value=freq_q2,
            min_val=5.595e9 ,
            max_val=5.605e9 ,
            unit='Hz 2pi'
        ),
        anhar=Qty(
            value=anhar_q2,
            min_val=-380e6 ,
            max_val=-120e6 ,
            unit='Hz 2pi'
        ),
        hilbert_dim=qubit_lvls,
        t1=Qty(
            value=t1_q2,
            min_val=1e-6,
            max_val=90e-6,
            unit='s'
        ),
        t2star=Qty(
            value=t2star_q2,
            min_val=10e-6,
            max_val=90e-6,
            unit='s'
        ),
        temp=Qty(
            value=qubit_temp,
            min_val=0.0,
            max_val=0.12,
            unit='K'
        )
    )

A static coupling between the two is realized in the following way. We
supply the type of coupling by selecting ``int_XX``
:math:`(b_1+b_1^\dagger)(b_2+b_2^\dagger)` from the hamiltonian library.
The “connected” property contains the list of qubit names to be coupled,
in this case “Q1” and “Q2”.

.. code-block:: python

    coupling_strength = 20e6
    q1q2 = chip.Coupling(
        name="Q1-Q2",
        desc="coupling",
        comment="Coupling qubit 1 to qubit 2",
        connected=["Q1", "Q2"],
        strength=Qty(
            value=coupling_strength,
            min_val=-1 * 1e3 ,
            max_val=200e6 ,
            unit='Hz 2pi'
        ),
        hamiltonian_func=hamiltonians.int_XX
    )

In the same spirit, we specify control Hamiltonians to drive the system.
Again “connected” connected tells us which qubit this drive acts on and
“name” will later be used to assign the correct control signal to this
drive line.

.. code-block:: python

    drive = chip.Drive(
        name="d1",
        desc="Drive 1",
        comment="Drive line 1 on qubit 1",
        connected=["Q1"],
        hamiltonian_func=hamiltonians.x_drive
    )
    drive2 = chip.Drive(
        name="d2",
        desc="Drive 2",
        comment="Drive line 2 on qubit 2",
        connected=["Q2"],
        hamiltonian_func=hamiltonians.x_drive
    )

SPAM errors
^^^^^^^^^^^

In experimental practice, the qubit state can be mis-classified during
read-out. We simulate this by constructing a *confusion matrix*,
containing the probabilities for one qubit state being mistaken for
another.

.. code-block:: python

    m00_q1 = 0.97 # Prop to read qubit 1 state 0 as 0
    m01_q1 = 0.04 # Prop to read qubit 1 state 0 as 1
    m00_q2 = 0.96 # Prop to read qubit 2 state 0 as 0
    m01_q2 = 0.05 # Prop to read qubit 2 state 0 as 1
    one_zeros = np.array([0] * qubit_lvls)
    zero_ones = np.array([1] * qubit_lvls)
    one_zeros[0] = 1
    zero_ones[0] = 0
    val1 = one_zeros * m00_q1 + zero_ones * m01_q1
    val2 = one_zeros * m00_q2 + zero_ones * m01_q2
    min_val = one_zeros * 0.8 + zero_ones * 0.0
    max_val = one_zeros * 1.0 + zero_ones * 0.2
    confusion_row1 = Qty(value=val1, min_val=min_val, max_val=max_val, unit="")
    confusion_row2 = Qty(value=val2, min_val=min_val, max_val=max_val, unit="")
    conf_matrix = tasks.ConfusionMatrix(Q1=confusion_row1, Q2=confusion_row2)

The following task creates an initial thermal state with given
temperature.

.. code-block:: python

    init_temp = 50e-3
    init_ground = tasks.InitialiseGround(
        init_temp=Qty(
            value=init_temp,
            min_val=-0.001,
            max_val=0.22,
            unit='K'
        )
    )

We collect the parts specified above in the Model.

.. code-block:: python

    model = Mdl(
        [q1, q2], # Individual, self-contained components
        [drive, drive2, q1q2],  # Interactions between components
        [conf_matrix, init_ground] # SPAM processing
    )

Further, we can decide between coherent or open-system dynamics using
set_lindbladian() and whether to eliminate the static coupling by going
to the dressed frame with set_dressed().

.. code:: ipython3

    model.set_lindbladian(False)
    model.set_dressed(True)

Control signals
^^^^^^^^^^^^^^^

With the system model taken care of, we now specify the control
electronics and signal chain. Complex shaped controls are often realized
by creating an envelope signal with an arbitrary waveform generator
(AWG) with limited bandwith and mixing it with a fast, stable local
oscillator (LO).

.. code-block:: python

    sim_res = 100e9 # Resolution for numerical simulation
    awg_res = 2e9 # Realistic, limited resolution of an AWG
    lo = devices.LO(name='lo', resolution=sim_res)
    awg = devices.AWG(name='awg', resolution=awg_res)
    mixer = devices.Mixer(name='mixer')

Waveform generators exhibit a rise time, the time it takes until the
target voltage is set. This has a smoothing effect on the resulting
pulse shape.

.. code-block:: python

    resp = devices.Response(
        name='resp',
        rise_time=Qty(
            value=0.3e-9,
            min_val=0.05e-9,
            max_val=0.6e-9,
            unit='s'
        ),
        resolution=sim_res
    )

In simulation, we translate between AWG resolution and simulation (or
“analog”) resolution by including an up-sampling device.

.. code-block:: python

    dig_to_an = devices.DigitalToAnalog(
        name="dac",
        resolution=sim_res
    )

Control electronics apply voltages to lines, whereas in a Hamiltonian we
usually write the control fields in energy or frequency units. In
practice, this conversion can be highly non-trivial if it involves
multiple stages of attenuation and for example the conversion of a line
voltage in an antenna to a dipole field coupling to the qubit. The
following device represents a simple, linear conversion factor.

.. code-block:: python

    v2hz = 1e9
    v_to_hz = devices.VoltsToHertz(
        name='v_to_hz',
        V_to_Hz=Qty(
            value=v2hz,
            min_val=0.9e9,
            max_val=1.1e9,
            unit='Hz/V'
        )
    )

The generator combines the parts of the signal generation and assigns a
signal chain to each control line.

.. code-block:: python

    generator = Gnr(
            devices={
                "LO": devices.LO(name='lo', resolution=sim_res, outputs=1),
                "AWG": devices.AWG(name='awg', resolution=awg_res, outputs=1),
                "DigitalToAnalog": devices.DigitalToAnalog(
                    name="dac",
                    resolution=sim_res,
                    inputs=1,
                    outputs=1
                ),
                "Response": devices.Response(
                    name='resp',
                    rise_time=Qty(
                        value=0.3e-9,
                        min_val=0.05e-9,
                        max_val=0.6e-9,
                        unit='s'
                    ),
                    resolution=sim_res,
                    inputs=1,
                    outputs=1
                ),
                "Mixer": devices.Mixer(name='mixer', inputs=2, outputs=1),
                "VoltsToHertz": devices.VoltsToHertz(
                    name='v_to_hz',
                    V_to_Hz=Qty(
                        value=1e9,
                        min_val=0.9e9,
                        max_val=1.1e9,
                        unit='Hz/V'
                    ),
                    inputs=1,
                    outputs=1
                )
            },
            chains= {
                "d1": {
                    "LO": [],
                    "AWG": [],
                    "DigitalToAnalog": ["AWG"],
                    "Response": ["DigitalToAnalog"],
                    "Mixer": ["LO", "Response"],
                    "VoltsToHertz": ["Mixer"]
                },
                "d2": {
                    "LO": [],
                    "AWG": [],
                    "DigitalToAnalog": ["AWG"],
                    "Response": ["DigitalToAnalog"],
                    "Mixer": ["LO", "Response"],
                    "VoltsToHertz": ["Mixer"]
                }
            }
        )

Gates-set and Parameter map
^^^^^^^^^^^^^^^^^^^^^^^^^^^

It remains to write down what kind of operations we want to perform on
the device. For a gate based quantum computing chip, we define a
gate-set.

We choose a gate time of 7ns and a Gaussian envelope shape with a list
of parameters.

.. code-block:: python

    t_final = 7e-9   # Time for single qubit gates
    sideband = 50e6
    gauss_params_single = {
        'amp': Qty(
            value=0.5,
            min_val=0.4,
            max_val=0.6,
            unit="V"
        ),
        't_final': Qty(
            value=t_final,
            min_val=0.5 * t_final,
            max_val=1.5 * t_final,
            unit="s"
        ),
        'sigma': Qty(
            value=t_final / 4,
            min_val=t_final / 8,
            max_val=t_final / 2,
            unit="s"
        ),
        'xy_angle': Qty(
            value=0.0,
            min_val=-0.5 * np.pi,
            max_val=2.5 * np.pi,
            unit='rad'
        ),
        'freq_offset': Qty(
            value=-sideband - 3e6 ,
            min_val=-56 * 1e6 ,
            max_val=-52 * 1e6 ,
            unit='Hz 2pi'
        ),
        'delta': Qty(
            value=-1,
            min_val=-5,
            max_val=3,
            unit=""
        )
    }

Here we take ``gaussian_nonorm()`` from the libraries as the function to
define the shape.

.. code-block:: python

    gauss_env_single = pulse.Envelope(
        name="gauss",
        desc="Gaussian comp for single-qubit gates",
        params=gauss_params_single,
        shape=envelopes.gaussian_nonorm
    )

We also define a gate that represents no driving.

.. code-block:: python

    nodrive_env = pulse.Envelope(
        name="no_drive",
        params={
            't_final': Qty(
                value=t_final,
                min_val=0.5 * t_final,
                max_val=1.5 * t_final,
                unit="s"
            )
        },
        shape=envelopes.no_drive
    )

We specify the drive tones with an offset from the qubit frequencies. As
is done in experiment, we will later adjust the resonance by modulating
the envelope function.

.. code-block:: python

    lo_freq_q1 = 5e9  + sideband
    carrier_parameters = {
        'freq': Qty(
            value=lo_freq_q1,
            min_val=4.5e9 ,
            max_val=6e9 ,
            unit='Hz 2pi'
        ),
        'framechange': Qty(
            value=0.0,
            min_val= -np.pi,
            max_val= 3 * np.pi,
            unit='rad'
        )
    }
    carr = pulse.Carrier(
        name="carrier",
        desc="Frequency of the local oscillator",
        params=carrier_parameters
    )

For the second qubit drive tone, we copy the first one and replace the
frequency. The deepcopy is to ensure that we don’t just create a pointer
to the first drive.

.. code-block:: python

    lo_freq_q2 = 5.6e9  + sideband
    carr_2 = copy.deepcopy(carr)
    carr_2.params['freq'].set_value(lo_freq_q2)

Instructions
^^^^^^^^^^^^

We define the gates we want to perform with a “name” that will identify
them later and “channels” relating to the control Hamiltonians and drive
lines we specified earlier. As a start we write down 90 degree rotations
in the positive :math:`x`-direction and identity gates for both qubits.
Then we add a carrier and envelope to each.

.. code-block:: python

    rx90p_q1 = gates.Instruction(
        name="rx90p", targets=[0], t_start=0.0, t_end=t_final, channels=["d1", "d2"]
    )
    rx90p_q2 = gates.Instruction(
        name="rx90p", targets=[1], t_start=0.0, t_end=t_final, channels=["d1", "d2"]
    )

    rx90p_q1.add_component(gauss_env_single, "d1")
    rx90p_q1.add_component(carr, "d1")


    rx90p_q2.add_component(copy.deepcopy(gauss_env_single), "d2")
    rx90p_q2.add_component(carr_2, "d2")


When later compiling gates into sequences, we have to take care of the
relative rotating frames of the qubits and local oscillators. We do this
by adding a phase after each gate that realigns the frames.

.. code-block:: python

    rx90p_q1.add_component(nodrive_env, "d2")
    rx90p_q1.add_component(copy.deepcopy(carr_2), "d2")
    rx90p_q1.comps["d2"]["carrier"].params["framechange"].set_value(
        (-sideband * t_final) * 2 * np.pi % (2 * np.pi)
    )

    rx90p_q2.add_component(nodrive_env, "d1")
    rx90p_q2.add_component(copy.deepcopy(carr), "d1")
    rx90p_q2.comps["d1"]["carrier"].params["framechange"].set_value(
        (-sideband * t_final) * 2 * np.pi % (2 * np.pi)
    )

The remainder of the gates-set can be derived from the RX90p gate by
shifting its phase by multiples of :math:`\pi/2`.

.. code-block:: python

    ry90p_q1 = copy.deepcopy(rx90p_q1)
    ry90p_q1.name = "ry90p"
    rx90m_q1 = copy.deepcopy(rx90p_q1)
    rx90m_q1.name = "rx90m"
    ry90m_q1 = copy.deepcopy(rx90p_q1)
    ry90m_q1.name = "ry90m"
    ry90p_q1.comps['d1']['gauss'].params['xy_angle'].set_value(0.5 * np.pi)
    rx90m_q1.comps['d1']['gauss'].params['xy_angle'].set_value(np.pi)
    ry90m_q1.comps['d1']['gauss'].params['xy_angle'].set_value(1.5 * np.pi)
    single_q_gates = [rx90p_q1, ry90p_q1, rx90m_q1, ry90m_q1]


    ry90p_q2 = copy.deepcopy(rx90p_q2)
    ry90p_q2.name = "ry90p"
    rx90m_q2 = copy.deepcopy(rx90p_q2)
    rx90m_q2.name = "rx90m"
    ry90m_q2 = copy.deepcopy(rx90p_q2)
    ry90m_q2.name = "ry90m"
    ry90p_q2.comps['d2']['gauss'].params['xy_angle'].set_value(0.5 * np.pi)
    rx90m_q2.comps['d2']['gauss'].params['xy_angle'].set_value(np.pi)
    ry90m_q2.comps['d2']['gauss'].params['xy_angle'].set_value(1.5 * np.pi)
    single_q_gates.extend([rx90p_q2, ry90p_q2, rx90m_q2, ry90m_q2])

With every component defined, we collect them in the parameter map, our
object that holds information and methods to manipulate and examine
model and control parameters.

.. code-block:: python

    parameter_map = PMap(instructions=all_1q_gates_comb, model=model, generator=generator)

The experiment
^^^^^^^^^^^^^^

Finally everything is collected in the experiment that provides the
functions to interact with the system.

.. code-block:: python

    exp = Exp(pmap=parameter_map)

Simulation
^^^^^^^^^^

With our experiment all set-up, we can perform simulations. We first
decide which basic gates to simulate, in this case only the 90 degree
rotation on one qubit and the identity.

.. code-block:: python

    exp.set_opt_gates(['RX90p:Id', 'Id:Id'])

The actual numerical simulation is done by calling ``exp.compute_propagators()``.
This is the most resource intensive part as it involves solving the
equations of motion for the system.

.. code-block:: python

    unitaries = exp.compute_propagators()


.. code-block:: python



Dynamics
~~~~~~~~

To investigate dynamics, we define the ground state as an initial state.

.. code:: ipython3

    psi_init = [[0] * 9]
    psi_init[0][0] = 1
    init_state = tf.transpose(tf.constant(psi_init, tf.complex128))

.. code:: ipython3

    init_state




.. parsed-literal::

    <tf.Tensor: shape=(9, 1), dtype=complex128, numpy=
    array([[1.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j],
           [0.+0.j]])>



Since we stored the process matrices, we can now relatively inexpensively
evaluate sequences. We start with just one gate

.. code:: ipython3

    barely_a_seq = ['rx90p[0]']

and plot system dynamics.

.. code-block:: python

    def plot_dynamics(exp, psi_init, seq, goal=-1):
            """
            Plotting code for time-resolved populations.

            Parameters
            ----------
            psi_init: tf.Tensor
                Initial state or density matrix.
            seq: list
                List of operations to apply to the initial state.
            goal: tf.float64
                Value of the goal function, if used.
            debug: boolean
                If true, return a matplotlib figure instead of saving.
            """
            model = exp.pmap.model
            dUs = exp.partial_propagators
            psi_t = psi_init.numpy()
            pop_t = exp.populations(psi_t, model.lindbladian)
            for gate in seq:
                for du in dUs[gate]:
                    psi_t = np.matmul(du.numpy(), psi_t)
                    pops = exp.populations(psi_t, model.lindbladian)
                    pop_t = np.append(pop_t, pops, axis=1)

            fig, axs = plt.subplots(1, 1)
            ts = exp.ts
            dt = ts[1] - ts[0]
            ts = np.linspace(0.0, dt*pop_t.shape[1], pop_t.shape[1])
            axs.plot(ts / 1e-9, pop_t.T)
            axs.grid(linestyle="--")
            axs.tick_params(
                direction="in", left=True, right=True, top=True, bottom=True
            )
            axs.set_xlabel('Time [ns]')
            axs.set_ylabel('Population')
            plt.legend(model.state_labels)
            pass



.. code-block:: python

    plot_dynamics(exp, init_state, barely_a_seq)



.. image:: dyn_singleX.png


We can see an ill-defined un-optimized gate. The labels indicate qubit
states in the product basis. Next we increase the number of repetitions
of the same gate.

.. code-block:: python

    barely_a_seq * 10




.. parsed-literal::

    ['rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]',
     'rx90p[0]']



.. code-block:: python

    plot_dynamics(exp, init_state, barely_a_seq * 5)



.. image:: dyn_5X.png


.. code-block:: python

    plot_dynamics(exp, init_state, barely_a_seq * 10)



.. image:: dyn_10X.png


Note that at this point, we only multiply already computed matrices. We
don’t need to solve the equations of motion again for new sequences.