Alexander Pitchford (agp1@aber.ac.uk)
Example to demonstrate using 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 initial ctrl pulses. Also re-uses objects in repeated runs with different total evolution times.
The (default) L-BFGS-B algorithm is used to optimise the pulse to minimise the fidelity error, which is equivalent maximising the fidelity to optimal value of 1.
The system in this example is two qubits in constant fields in x, y and z with 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:
The user can experiment with the timeslicing, by means of changing the timeslots durations. Different initial (starting) pulse types can be tried. The initial and final pulses are displayed in a plot
This example assumes that the example-control-pulseoptim-Hadamard has already been tried, and hence explanations in that notebook are not repeated here.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import datetime
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'
Note here that there are two controls acting on each qubit.
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)
Multiple total evolution times will be tried. Using this approach, the minimum evolution time required to achieve the target fidelity could be determined (iteratively).
Note that the timeslot duration dt is fixed, and so the number of timeslots depends on the evo_time
# Duration of each timeslot
dt = 0.05
# List of evolution times to try
evo_times = [1, 3, 6]
n_evo_times = len(evo_times)
evo_time = evo_times[0]
n_ts = int(float(evo_time) / dt)
#Empty list that will hold the results for each evolution time
results = list()
# Fidelity error target
fid_err_targ = 1e-5
# Maximum iterations for the optisation algorithm
max_iter = 200
# Maximum (elapsed) time allowed in seconds
max_wall_time = 120
# Minimum gradient (sum of gradients squared)
# as this tends to 0 -> local minima has been found
min_grad = 1e-20
Here the linear initial pulse type is used, simply because it results in smooth final pulses
# pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|
p_type = 'LIN'
#Set to None to suppress output files
f_ext = "{}_n_ts{}_ptype{}.txt".format(example_name, n_ts, p_type)
Here is the main difference between this and the Hadamard example. In this case we use a different pulseoptim function that just creates the objects that can be used to set the physics and configure the optimisation algorithm. This gives greater flexibility (shown here by seting different initial pulse parameters for each control) and is also more efficient when running multiple optimisations on the same system.
optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ, n_ts, evo_time,
amp_lbound=-5.0, amp_ubound=5.0,
fid_err_targ=fid_err_targ, min_grad=min_grad,
max_iter=max_iter, max_wall_time=max_wall_time,
optim_method='fmin_l_bfgs_b',
method_params={'max_metric_corr':20, 'accuracy_factor':1e8},
dyn_type='UNIT',
fid_params={'phase_option':'PSU'},
init_pulse_type=p_type,
log_level=log_level, gen_stats=True)
# **** get handles to the other objects ****
optim.test_out_files = 0
dyn = optim.dynamics
dyn.test_out_files = 0
p_gen = optim.pulse_generator
Here a loop is used to perform the optimisation for each of the evo_times given in the list above. The first optimisation is completed using the timeslot parameters passed when the optimisation objects are created. For the subsequent runs, the Dynamics object 'dyn' is used to set the timeslot parameters before the initial pulses are generated and optimisation is completed. Note that using this method, the dyn.initialize_controls method must be called with an array of the initial amplitudes before the optim.run_optimization method is called.
for i in range(n_evo_times):
# Generate the tau (duration) and time (cumulative) arrays
# so that it can be used to create the pulse generator
# with matching timeslots
dyn.init_timeslots()
if i > 0:
# Create a new pulse generator for the new dynamics
p_gen = pulsegen.create_pulse_gen(p_type, dyn)
#Generate different initial pulses for each of the controls
init_amps = np.zeros([n_ts, n_ctrls])
if (p_gen.periodic):
phase_diff = np.pi / n_ctrls
for j in range(n_ctrls):
init_amps[:, j] = p_gen.gen_pulse(start_phase=phase_diff*j)
elif (isinstance(p_gen, pulsegen.PulseGenLinear)):
for j in range(n_ctrls):
p_gen.scaling = float(j) - float(n_ctrls - 1)/2
init_amps[:, j] = p_gen.gen_pulse()
elif (isinstance(p_gen, pulsegen.PulseGenZero)):
for j in range(n_ctrls):
p_gen.offset = sf = float(j) - float(n_ctrls - 1)/2
init_amps[:, j] = p_gen.gen_pulse()
else:
# Should be random pulse
for j in range(n_ctrls):
init_amps[:, j] = p_gen.gen_pulse()
dyn.initialize_controls(init_amps)
# 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("\n+++++++++++++++++++++++++++++++++++")
print("Starting pulse optimisation for T={}".format(evo_time))
print("+++++++++++++++++++++++++++++++++++\n")
result = optim.run_optimization()
results.append(result)
# 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)
# Report the results
result.stats.report()
print("Final evolution\n{}\n".format(result.evo_full_final))
print("********* Summary *****************")
print("Final fidelity error {}".format(result.fid_err))
print("Final gradient normal {}".format(result.grad_norm_final))
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)))
if i+1 < len(evo_times):
# reconfigure the dynamics for the next evo time
evo_time = evo_times[i+1]
n_ts = int(float(evo_time) / dt)
dyn.tau = None
dyn.evo_time = evo_time
dyn.num_tslots = n_ts
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 INFO:qutip.control.optimizer:Optimising pulse(s) using GRAPE with 'fmin_l_bfgs_b' method
Initial amplitudes output to file: ctrl_amps_initial_QFT_n_ts20_ptypeLIN.txt *********************************** +++++++++++++++++++++++++++++++++++ Starting pulse optimisation for T=1 +++++++++++++++++++++++++++++++++++
INFO:qutip.control.dynamics:Setting memory optimisations for level 0 INFO:qutip.control.dynamics:Using operator data type <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 INFO:qutip.control.optimizer:Optimising pulse(s) using GRAPE with 'fmin_l_bfgs_b' method
Final amplitudes output to file: ctrl_amps_final_QFT_n_ts20_ptypeLIN.txt ------------------------------------ ---- Control optimisation stats ---- **** Timings (HH:MM:SS.US) **** Total wall time elapsed during optimisation: 0:00:00.277534 Wall time computing Hamiltonians: 0:00:00.015000 (5.40%) Wall time computing propagators: 0:00:00.191946 (69.16%) Wall time computing forward propagation: 0:00:00.002414 (0.87%) Wall time computing onward propagation: 0:00:00.002215 (0.80%) Wall time computing gradient: 0:00:00.044586 (16.07%) **** Iterations and function calls **** Number of iterations: 60 Number of fidelity function calls: 64 Number of times fidelity is computed: 64 Number of gradient function calls: 64 Number of times gradients are computed: 64 Number of times timeslot evolution is recomputed: 64 **** Control amplitudes **** Number of control amplitude updates: 63 Mean number of updates per iteration: 1.05 Number of timeslot values changed: 701 Mean number of timeslot changes per update: 11.126984126984127 Number of amplitude values changed: 1343 Mean number of amplitude changes per update: 21.317460317460316 ------------------------------------ Final evolution Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False Qobj data = [[-0.16998421+0.48769402j 0.00937949+0.31030265j -0.78519291+0.0495552j -0.06609584+0.11632679j] [ 0.01430120+0.29484354j 0.12448164-0.48737042j -0.07997799-0.24264738j 0.76748036+0.07440979j] [-0.79272858+0.02269328j -0.03496912-0.24748526j 0.08419464+0.50446146j 0.04342906-0.21245808j] [-0.03607690+0.12604581j 0.76668361-0.01797591j 0.00951480-0.23255432j -0.20301010-0.54708215j]] ********* Summary ***************** Final fidelity error 0.3436799944999004 Final gradient normal 0.015892258406309603 Terminated due to function converged Number of iterations 60 Completed in 0:00:00.277534 HH:MM:SS.US Initial amplitudes output to file: ctrl_amps_initial_QFT_n_ts20_ptypeLIN.txt *********************************** +++++++++++++++++++++++++++++++++++ Starting pulse optimisation for T=3 +++++++++++++++++++++++++++++++++++
INFO:qutip.control.dynamics:Setting memory optimisations for level 0 INFO:qutip.control.dynamics:Using operator data type <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 INFO:qutip.control.optimizer:Optimising pulse(s) using GRAPE with 'fmin_l_bfgs_b' method
Final amplitudes output to file: ctrl_amps_final_QFT_n_ts20_ptypeLIN.txt ------------------------------------ ---- Control optimisation stats ---- **** Timings (HH:MM:SS.US) **** Total wall time elapsed during optimisation: 0:00:00.390896 Wall time computing Hamiltonians: 0:00:00.023207 (5.94%) Wall time computing propagators: 0:00:00.279448 (71.49%) Wall time computing forward propagation: 0:00:00.003514 (0.90%) Wall time computing onward propagation: 0:00:00.003310 (0.85%) Wall time computing gradient: 0:00:00.068349 (17.49%) **** Iterations and function calls **** Number of iterations: 29 Number of fidelity function calls: 36 Number of times fidelity is computed: 36 Number of gradient function calls: 35 Number of times gradients are computed: 35 Number of times timeslot evolution is recomputed: 36 **** Control amplitudes **** Number of control amplitude updates: 35 Mean number of updates per iteration: 1.206896551724138 Number of timeslot values changed: 2100 Mean number of timeslot changes per update: 60.0 Number of amplitude values changed: 8390 Mean number of amplitude changes per update: 239.71428571428572 ------------------------------------ Final evolution Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False Qobj data = [[-0.18922387+0.46120317j -0.19352353+0.46212201j -0.19278561+0.46360969j -0.18771204+0.46167283j] [-0.18926712+0.46237894j -0.46052942-0.19175575j 0.19104394-0.46161315j 0.46401084+0.19140406j] [-0.19133929+0.46005139j 0.19127929-0.46261849j -0.19190677+0.46291531j 0.19214148-0.46162888j] [-0.19289196+0.46519277j 0.46053463+0.1934903j 0.19031547-0.45932179j -0.46156574-0.19135913j]] ********* Summary ***************** Final fidelity error 9.093450580754947e-06 Final gradient normal 0.00021946915669970356 Terminated due to Goal achieved Number of iterations 29 Completed in 0:00:00.390896 HH:MM:SS.US Initial amplitudes output to file: ctrl_amps_initial_QFT_n_ts20_ptypeLIN.txt *********************************** +++++++++++++++++++++++++++++++++++ Starting pulse optimisation for T=6 +++++++++++++++++++++++++++++++++++ Final amplitudes output to file: ctrl_amps_final_QFT_n_ts20_ptypeLIN.txt ------------------------------------ ---- Control optimisation stats ---- **** Timings (HH:MM:SS.US) **** Total wall time elapsed during optimisation: 0:00:00.650126 Wall time computing Hamiltonians: 0:00:00.041716 (6.42%) Wall time computing propagators: 0:00:00.468142 (72.01%) Wall time computing forward propagation: 0:00:00.005879 (0.90%) Wall time computing onward propagation: 0:00:00.005721 (0.88%) Wall time computing gradient: 0:00:00.115748 (17.80%) **** Iterations and function calls **** Number of iterations: 24 Number of fidelity function calls: 30 Number of times fidelity is computed: 30 Number of gradient function calls: 29 Number of times gradients are computed: 29 Number of times timeslot evolution is recomputed: 30 **** Control amplitudes **** Number of control amplitude updates: 29 Mean number of updates per iteration: 1.2083333333333333 Number of timeslot values changed: 3480 Mean number of timeslot changes per update: 120.0 Number of amplitude values changed: 13920 Mean number of amplitude changes per update: 480.0 ------------------------------------ Final evolution Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False Qobj data = [[ 0.46207580+0.19122715j 0.46151636+0.19126417j 0.46241052+0.19094127j 0.46213112+0.19102767j] [ 0.46199149+0.19151454j -0.19197958+0.46204899j -0.46180318-0.19115039j 0.19122368-0.46170712j] [ 0.46165031+0.19097801j -0.46193761-0.1917004j 0.46197027+0.19137927j -0.46218878-0.19133733j] [ 0.46216726+0.19134259j 0.19167781-0.46173526j -0.46181095-0.19132534j -0.19139753+0.46188943j]] ********* Summary ***************** Final fidelity error 2.4336489157228414e-07 Final gradient normal 0.0006194193096196838 Terminated due to Goal achieved Number of iterations 24 Completed in 0:00:00.650126 HH:MM:SS.US
fig1 = plt.figure(figsize=(12,8))
for i in range(n_evo_times):
#Initial amps
ax1 = fig1.add_subplot(2, n_evo_times, i+1)
ax1.set_title("Init amps T={}".format(evo_times[i]))
# ax1.set_xlabel("Time")
ax1.get_xaxis().set_visible(False)
if i == 0:
ax1.set_ylabel("Control amplitude")
for j in range(n_ctrls):
ax1.step(results[i].time,
np.hstack((results[i].initial_amps[:, j],
results[i].initial_amps[-1, j])),
where='post')
ax2 = fig1.add_subplot(2, n_evo_times, i+n_evo_times+1)
ax2.set_title("Final amps T={}".format(evo_times[i]))
ax2.set_xlabel("Time")
#Optimised amps
if i == 0:
ax2.set_ylabel("Control amplitude")
for j in range(n_ctrls):
ax2.step(results[i].time,
np.hstack((results[i].final_amps[:, j],
results[i].final_amps[-1, j])),
where='post')
plt.tight_layout()
plt.show()
from qutip.ipynbtools import version_table
version_table()
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 17:13:52 2017 BST |