Notebook

Calculation of control fields for QFT gate on two qubits using the CRAB algorithm¶

Alexander Pitchford (agp1@aber.ac.uk)

Example to demonstrate using the CRAB [1][2] algorithm in the control library to determine control pulses using the ctrlpulseoptim.create_pulse_optimizer function to generate an Optimizer object, through which the configuration can be manipulated before running the optmisation algorithm. In this case it is demonstrated by modifying the CRAB pulse parameters to show how pulse constraints for controls can be applied. The system in this example is two qubits in constant fields in x, y and z with a variable independant controls fields in x and y acting on each qubit The target evolution is the QFT gate. The user can experiment with the different: phase options - phase_option = SU or PSU propagtor computer type prop_type = DIAG or FRECHET fidelity measures - fid_type = UNIT or TRACEDIFF The user can experiment with the timeslicing, by means of changing the number of timeslots and/or total time for the evolution. Different guess and ramping pulse parameters can be tried. The initial and final pulses are displayed in a plot

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import datetime

In [2]:

from qutip import Qobj, identity, sigmax, sigmay, sigmaz, tensor
from qutip.qip.algorithms import qft
import qutip.logging_utils as logging
logger = logging.get_logger()
#Set this to None or logging.WARN for 'quiet' execution
log_level = logging.INFO
#QuTiP control modules
import qutip.control.pulseoptim as cpo
import qutip.control.pulsegen as pulsegen

example_name = 'QFT'

Defining the physics¶

In [3]:

Sx = sigmax()
Sy = sigmay()
Sz = sigmaz()
Si = 0.5*identity(2)

# Drift Hamiltonian
H_d = 0.5*(tensor(Sx, Sx) + tensor(Sy, Sy) + tensor(Sz, Sz))
# The (four) control Hamiltonians
H_c = [tensor(Sx, Si), tensor(Sy, Si), tensor(Si, Sx), tensor(Si, Sy)]
n_ctrls = len(H_c)
# start point for the gate evolution
U_0 = identity(4)
# Target for the gate evolution - Quantum Fourier Transform gate
U_targ = qft.qft(2)

Defining the time evolution parameters¶

In [4]:

# Number of time slots
n_ts = 200
# Time allowed for the evolution
evo_time = 10

Set the conditions which will cause the pulse optimisation to terminate¶

In [5]:

# Fidelity error target
fid_err_targ = 1e-3
# Maximum iterations for the optisation algorithm
max_iter = 20000
# Maximum (elapsed) time allowed in seconds
max_wall_time = 300

Give an extension for output files¶

In [6]:

#Set to None to suppress output files
f_ext = "{}_n_ts{}.txt".format(example_name, n_ts)

Create the optimiser objects¶

In [7]:

optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ, n_ts, evo_time, 
                fid_err_targ=fid_err_targ, 
                max_iter=max_iter, max_wall_time=max_wall_time,
                alg='CRAB', 
                dyn_type='UNIT', 
                prop_type='DIAG', 
                fid_type='UNIT', fid_params={'phase_option':'PSU'}, 
                log_level=log_level, gen_stats=True)
                

Configure the pulses for each of the controls¶

In [8]:

dyn = optim.dynamics

# Control 1
crab_pgen = optim.pulse_generator[0]
# Start from a ramped pulse
guess_pgen = pulsegen.create_pulse_gen('LIN', dyn=dyn, 
                                           pulse_params={'scaling':3.0})
crab_pgen.guess_pulse = guess_pgen.gen_pulse()
crab_pgen.scaling = 0.0
# Add some higher frequency components
crab_pgen.num_coeffs = 5

# Control 2
crab_pgen = optim.pulse_generator[1]
# Apply a ramping pulse that will force the start and end to zero
ramp_pgen = pulsegen.create_pulse_gen('GAUSSIAN_EDGE', dyn=dyn, 
                                    pulse_params={'decay_time':evo_time/50.0})
crab_pgen.ramping_pulse = ramp_pgen.gen_pulse()

# Control 3
crab_pgen = optim.pulse_generator[2]
# Add bounds
crab_pgen.scaling = 0.5
crab_pgen.lbound = -2.0
crab_pgen.ubound = 2.0


# Control 4
crab_pgen = optim.pulse_generator[3]
# Start from a triangular pulse with small signal
guess_pgen = pulsegen.PulseGenTriangle(dyn=dyn)
guess_pgen.num_waves = 1
guess_pgen.scaling = 2.0
guess_pgen.offset = 2.0
crab_pgen.guess_pulse = guess_pgen.gen_pulse()
crab_pgen.scaling = 0.1

init_amps = np.zeros([n_ts, n_ctrls])
for j in range(dyn.num_ctrls):
    pgen = optim.pulse_generator[j]
    pgen.init_pulse()
    init_amps[:, j] = pgen.gen_pulse()

dyn.initialize_controls(init_amps)

INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system
INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system
INFO:qutip.control.pulsegen:The number of CRAB coefficients per basis function has been estimated as 3, which means a total of 6 optimisation variables for this pulse. Based on the dimension (4) of the system
INFO:qutip.control.dynamics:Setting memory optimisations for level 0
INFO:qutip.control.dynamics:Internal operator data type choosen to be <class 'numpy.ndarray'>
INFO:qutip.control.dynamics:phased dynamics generator caching True
INFO:qutip.control.dynamics:propagator gradient caching True
INFO:qutip.control.dynamics:eigenvector adjoint caching True
INFO:qutip.control.dynamics:use sparse eigen decomp False

Run the pulse optimisation¶

In [9]:

# Save initial amplitudes to a text file
if f_ext is not None:
    pulsefile = "ctrl_amps_initial_" + f_ext
    dyn.save_amps(pulsefile)
    print("Initial amplitudes output to file: " + pulsefile)

print("***********************************")
print("Starting pulse optimisation")
result = optim.run_optimization()

# Save final amplitudes to a text file
if f_ext is not None:
    pulsefile = "ctrl_amps_final_" + f_ext
    dyn.save_amps(pulsefile)
    print("Final amplitudes output to file: " + pulsefile)

Initial amplitudes output to file: ctrl_amps_initial_QFT_n_ts200.txt
***********************************

INFO:qutip.control.optimizer:Optimising pulse(s) using CRAB with 'fmin' (Nelder-Mead) method

Starting pulse optimisation
Final amplitudes output to file: ctrl_amps_final_QFT_n_ts200.txt

Report the results¶

In [10]:

result.stats.report()
print("Final evolution\n{}\n".format(result.evo_full_final))
print("********* Summary *****************")
print("Initial fidelity error {}".format(result.initial_fid_err))
print("Final fidelity error {}".format(result.fid_err))
print("Terminated due to {}".format(result.termination_reason))
print("Number of iterations {}".format(result.num_iter))
print("Completed in {} HH:MM:SS.US".format(
        datetime.timedelta(seconds=result.wall_time)))

------------------------------------
---- Control optimisation stats ----
**** Timings (HH:MM:SS.US) ****
Total wall time elapsed during optimisation: 0:05:00.003775
Wall time computing Hamiltonians: 0:00:22.547817 (7.52%)
Wall time computing propagators: 0:04:24.789280 (88.26%)
Wall time computing forward propagation: 0:00:03.857988 (1.29%)
Wall time computing onward propagation: 0:00:03.634444 (1.21%)
Wall time computing gradient: 0:00:00 (0.00%)

**** Iterations and function calls ****
Number of iterations: 8610
Number of fidelity function calls: 10459
Number of times fidelity is computed: 10459
Number of gradient function calls: 0
Number of times gradients are computed: 0
Number of times timeslot evolution is recomputed: 10459

**** Control amplitudes ****
Number of control amplitude updates: 10458
Mean number of updates per iteration: 1.2146341463414634
Number of timeslot values changed: 2091599
Mean number of timeslot changes per update: 199.99990437942245
Number of amplitude values changed: 8350196
Mean number of amplitude changes per update: 798.4505641614076
------------------------------------
Final evolution
Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False
Qobj data =
[[-0.47890815-0.25124383j -0.47862720-0.14293885j -0.47089104-0.19654053j
  -0.39814103-0.19780087j]
 [-0.48524763-0.18754841j  0.19902638-0.45262542j  0.45666538+0.17490099j
  -0.22150294+0.44348832j]
 [-0.44119050-0.18412439j  0.49542987+0.17619236j -0.42828736-0.17952947j
   0.50360656+0.16023168j]
 [-0.41755379-0.18434167j -0.18861574+0.44037804j  0.50645054+0.16836502j
   0.23840544-0.4695553j ]]

********* Summary *****************
Initial fidelity error 0.9483036893449602
Final fidelity error 0.0032564547889728512
Terminated due to Max wall time exceeded
Number of iterations 8611
Completed in 0:05:00.003775 HH:MM:SS.US

Plot the initial and final amplitudes¶

In [12]:

fig1 = plt.figure()
ax1 = fig1.add_subplot(2, 1, 1)
ax1.set_title("Initial Control amps")
ax1.set_xlabel("Time")
ax1.set_ylabel("Control amplitude")
for j in range(n_ctrls):
    ax1.step(result.time, 
             np.hstack((result.initial_amps[:, j], result.initial_amps[-1, j])), 
             where='post')
ax2 = fig1.add_subplot(2, 1, 2)
ax2.set_title("Optimised Control Amplitudes")
ax2.set_xlabel("Time")
ax2.set_ylabel("Control amplitude")
for j in range(n_ctrls):
    ax2.step(result.time, 
             np.hstack((result.final_amps[:, j], result.final_amps[-1, j])), 
             where='post', label='u{}'.format(j))
ax2.legend(loc=8, ncol=n_ctrls)
plt.tight_layout()
plt.show()

Versions¶

In [13]:

from qutip.ipynbtools import version_table

version_table()

Out[13]:

Software	Version
QuTiP	4.1.0
Numpy	1.11.3
SciPy	0.18.1
matplotlib	2.0.0
Cython	0.25.2
Number of CPUs	4
BLAS Info	INTEL MKL
IPython	5.1.0
Python	3.6.0 \|Anaconda 4.3.1 (64-bit)\| (default, Dec 23 2016, 12:22:00) [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
OS	posix [linux]
Fri Jul 14 16:41:51 2017 BST

References: 3. Doria, P., Calarco, T. & Montangero, S. Optimal Control Technique for Many-Body Quantum Dynamics. Phys. Rev. Lett. 106, 1–4 (2011). 4. Caneva, T., Calarco, T. & Montangero, S. Chopped random-basis quantum optimization. Phys. Rev. A - At. Mol. Opt. Phys. 84, (2011).