# -*- coding: utf-8 -*-
# This file is part of QuTiP: Quantum Toolbox in Python.
#
# Copyright (c) 2014 and later, Alexander J G Pitchford
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
# of its contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
# @author: Alexander Pitchford
# @email1: agp1@aber.ac.uk
# @email2: alex.pitchford@gmail.com
# @organization: Aberystwyth University
# @supervisor: Daniel Burgarth
"""
Pulse generator - Generate pulses for the timeslots
Each class defines a gen_pulse function that produces a float array of
size num_tslots. Each class produces a differ type of pulse.
See the class and gen_pulse function descriptions for details
"""
import numpy as np
import qutip.logging_utils as logging
logger = logging.get_logger()
import qutip.control.dynamics as dynamics
import qutip.control.errors as errors
[docs]def create_pulse_gen(pulse_type='RND', dyn=None, pulse_params=None):
"""
Create and return a pulse generator object matching the given type.
The pulse generators each produce a different type of pulse,
see the gen_pulse function description for details.
These are the random pulse options:
RND - Independent random value in each timeslot
RNDFOURIER - Fourier series with random coefficients
RNDWAVES - Summation of random waves
RNDWALK1 - Random change in amplitude each timeslot
RNDWALK2 - Random change in amp gradient each timeslot
These are the other non-periodic options:
LIN - Linear, i.e. contant gradient over the time
ZERO - special case of the LIN pulse, where the gradient is 0
These are the periodic options
SINE - Sine wave
SQUARE - Square wave
SAW - Saw tooth wave
TRIANGLE - Triangular wave
If a Dynamics object is passed in then this is used in instantiate
the PulseGen, meaning that some timeslot and amplitude properties
are copied over.
"""
if pulse_type == 'RND':
return PulseGenRandom(dyn, params=pulse_params)
if pulse_type == 'RNDFOURIER':
return PulseGenRndFourier(dyn, params=pulse_params)
if pulse_type == 'RNDWAVES':
return PulseGenRndWaves(dyn, params=pulse_params)
if pulse_type == 'RNDWALK1':
return PulseGenRndWalk1(dyn, params=pulse_params)
if pulse_type == 'RNDWALK2':
return PulseGenRndWalk2(dyn, params=pulse_params)
elif pulse_type == 'LIN':
return PulseGenLinear(dyn, params=pulse_params)
elif pulse_type == 'ZERO':
return PulseGenZero(dyn, params=pulse_params)
elif pulse_type == 'SINE':
return PulseGenSine(dyn, params=pulse_params)
elif pulse_type == 'SQUARE':
return PulseGenSquare(dyn, params=pulse_params)
elif pulse_type == 'SAW':
return PulseGenSaw(dyn, params=pulse_params)
elif pulse_type == 'TRIANGLE':
return PulseGenTriangle(dyn, params=pulse_params)
elif pulse_type == 'GAUSSIAN':
return PulseGenGaussian(dyn, params=pulse_params)
elif pulse_type == 'CRAB_FOURIER':
return PulseGenCrabFourier(dyn, params=pulse_params)
elif pulse_type == 'GAUSSIAN_EDGE':
return PulseGenGaussianEdge(dyn, params=pulse_params)
else:
raise ValueError("No option for pulse_type '{}'".format(pulse_type))
[docs]class PulseGen(object):
"""
Pulse generator
Base class for all Pulse generators
The object can optionally be instantiated with a Dynamics object,
in which case the timeslots and amplitude scaling and offset
are copied from that.
Otherwise the class can be used independently by setting:
tau (array of timeslot durations)
or
num_tslots and pulse_time for equally spaced timeslots
Attributes
----------
num_tslots : integer
Number of timeslots, aka timeslices
(copied from Dynamics if given)
pulse_time : float
total duration of the pulse
(copied from Dynamics.evo_time if given)
scaling : float
linear scaling applied to the pulse
(copied from Dynamics.initial_ctrl_scaling if given)
offset : float
linear offset applied to the pulse
(copied from Dynamics.initial_ctrl_offset if given)
tau : array[num_tslots] of float
Duration of each timeslot
(copied from Dynamics if given)
lbound : float
Lower boundary for the pulse amplitudes
Note that the scaling and offset attributes can be used to fully
bound the pulse for all generators except some of the random ones
This bound (if set) may result in additional shifting / scaling
Default is -Inf
ubound : float
Upper boundary for the pulse amplitudes
Note that the scaling and offset attributes can be used to fully
bound the pulse for all generators except some of the random ones
This bound (if set) may result in additional shifting / scaling
Default is Inf
periodic : boolean
True if the pulse generator produces periodic pulses
random : boolean
True if the pulse generator produces random pulses
log_level : integer
level of messaging output from the logger.
Options are attributes of qutip.logging_utils,
in decreasing levels of messaging, are:
DEBUG_INTENSE, DEBUG_VERBOSE, DEBUG, INFO, WARN, ERROR, CRITICAL
Anything WARN or above is effectively 'quiet' execution,
assuming everything runs as expected.
The default NOTSET implies that the level will be taken from
the QuTiP settings file, which by default is WARN
"""
def __init__(self, dyn=None, params=None):
self.parent = dyn
self.params = params
self.reset()
[docs] def reset(self):
"""
reset attributes to default values
"""
if isinstance(self.parent, dynamics.Dynamics):
dyn = self.parent
self.num_tslots = dyn.num_tslots
self.pulse_time = dyn.evo_time
self.scaling = dyn.initial_ctrl_scaling
self.offset = dyn.initial_ctrl_offset
self.tau = dyn.tau
self.log_level = dyn.log_level
else:
self.num_tslots = 100
self.pulse_time = 1.0
self.scaling = 1.0
self.tau = None
self.offset = 0.0
self._uses_time = False
self.time = None
self._pulse_initialised = False
self.periodic = False
self.random = False
self.lbound = None
self.ubound = None
self.ramping_pulse = None
self.apply_params()
[docs] def apply_params(self, params=None):
"""
Set object attributes based on the dictionary (if any) passed in the
instantiation, or passed as a parameter
This is called during the instantiation automatically.
The key value pairs are the attribute name and value
"""
if not params:
params = self.params
if isinstance(params, dict):
self.params = params
for key in params:
setattr(self, key, params[key])
@property
def log_level(self):
return logger.level
@log_level.setter
def log_level(self, lvl):
"""
Set the log_level attribute and set the level of the logger
that is call logger.setLevel(lvl)
"""
logger.setLevel(lvl)
[docs] def gen_pulse(self):
"""
returns the pulse as an array of vales for each timeslot
Must be implemented by subclass
"""
# must be implemented by subclass
raise errors.UsageError(
"No method defined for generating a pulse. "
" Suspect base class was used where sub class should have been")
[docs] def init_pulse(self):
"""
Initialise the pulse parameters
"""
if self.tau is None:
self.tau = np.ones(self.num_tslots, dtype='f') * \
self.pulse_time/self.num_tslots
if self._uses_time:
self.time = np.zeros(self.num_tslots, dtype=float)
for k in range(self.num_tslots-1):
self.time[k+1] = self.time[k] + self.tau[k]
self._pulse_initialised = True
if not self.lbound is None:
if np.isinf(self.lbound):
self.lbound = None
if not self.ubound is None:
if np.isinf(self.ubound):
self.ubound = None
if not self.ubound is None and not self.lbound is None:
if self.ubound < self.lbound:
raise ValueError("ubound cannot be less the lbound")
def _apply_bounds_and_offset(self, pulse):
"""
Ensure that the randomly generated pulse fits within the bounds
(after applying the offset)
Assumes that pulses passed are centered around zero (on average)
"""
if self.lbound is None and self.ubound is None:
return pulse + self.offset
max_amp = max(pulse)
min_amp = min(pulse)
if ((self.ubound is None or max_amp + self.offset <= self.ubound) and
(self.lbound is None or min_amp + self.offset >= self.lbound)):
return pulse + self.offset
# Some shifting / scaling is required.
if self.ubound is None or self.lbound is None:
# One of the bounds is inf, so just shift the pulse
if self.lbound is None:
# max_amp + offset must exceed the ubound
return pulse + self.ubound - max_amp
else:
# min_amp + offset must exceed the lbound
return pulse + self.lbound - min_amp
else:
bound_range = self.ubound - self.lbound
amp_range = max_amp - min_amp
if max_amp - min_amp > bound_range:
# pulse range is too high, it must be scaled
pulse = pulse * bound_range / amp_range
# otherwise the pulse should fit anyway
return pulse + self.lbound - min(pulse)
def _apply_ramping_pulse(self, pulse, ramping_pulse=None):
if ramping_pulse is None:
ramping_pulse = self.ramping_pulse
if ramping_pulse is not None:
pulse = pulse*ramping_pulse
return pulse
[docs]class PulseGenZero(PulseGen):
"""
Generates a flat pulse
"""
[docs] def gen_pulse(self):
"""
Generate a pulse with the same value in every timeslot.
The value will be zero, unless the offset is not zero,
in which case it will be the offset
"""
pulse = np.zeros(self.num_tslots)
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenRandom(PulseGen):
"""
Generates random pulses as simply random values for each timeslot
"""
def reset(self):
PulseGen.reset(self)
self.random = True
self.apply_params()
[docs] def gen_pulse(self):
"""
Generate a pulse of random values between 1 and -1
Values are scaled using the scaling property
and shifted using the offset property
Returns the pulse as an array of vales for each timeslot
"""
pulse = (2*np.random.random(self.num_tslots) - 1) * self.scaling
return self._apply_bounds_and_offset(pulse)
class PulseGenRndFourier(PulseGen):
"""
Generates pulses by summing sine waves as a Fourier series
with random coefficients
Attributes
----------
scaling : float
The pulses should fit approximately within -/+scaling
(before the offset is applied)
as it is used to set a maximum for each component wave
Use bounds to be sure
(copied from Dynamics.initial_ctrl_scaling if given)
min_wavelen : float
Minimum wavelength of any component wave
Set by default to 1/10th of the pulse time
"""
def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.random = True
self._uses_time = True
try:
self.min_wavelen = self.pulse_time / 10.0
except:
self.min_wavelen = 0.1
self.apply_params()
def gen_pulse(self, min_wavelen=None):
"""
Generate a random pulse based on a Fourier series with a minimum
wavelength
"""
if min_wavelen is not None:
self.min_wavelen = min_wavelen
min_wavelen = self.min_wavelen
if min_wavelen > self.pulse_time:
raise ValueError("Minimum wavelength cannot be greater than "
"the pulse time")
if not self._pulse_initialised:
self.init_pulse()
# use some phase to avoid the first pulse being always 0
sum_wave = np.zeros(self.tau.shape)
wavelen = 2.0*self.pulse_time
t = self.time
wl = []
while wavelen > min_wavelen:
wl.append(wavelen)
wavelen = wavelen/2.0
num_comp_waves = len(wl)
amp_scale = np.sqrt(8)*self.scaling / float(num_comp_waves)
for wavelen in wl:
amp = amp_scale*(np.random.rand()*2 - 1)
phase_off = np.random.rand()*np.pi/2.0
curr_wave = amp*np.sin(2*np.pi*t/wavelen + phase_off)
sum_wave += curr_wave
return self._apply_bounds_and_offset(sum_wave)
class PulseGenRndWaves(PulseGen):
"""
Generates pulses by summing sine waves with random frequencies
amplitudes and phase offset
Attributes
----------
scaling : float
The pulses should fit approximately within -/+scaling
(before the offset is applied)
as it is used to set a maximum for each component wave
Use bounds to be sure
(copied from Dynamics.initial_ctrl_scaling if given)
num_comp_waves : integer
Number of component waves. That is the number of waves that
are summed to make the pulse signal
Set to 20 by default.
min_wavelen : float
Minimum wavelength of any component wave
Set by default to 1/10th of the pulse time
max_wavelen : float
Maximum wavelength of any component wave
Set by default to twice the pulse time
"""
def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.random = True
self._uses_time = True
self.num_comp_waves = 20
try:
self.min_wavelen = self.pulse_time / 10.0
except:
self.min_wavelen = 0.1
try:
self.max_wavelen = 2*self.pulse_time
except:
self.max_wavelen = 10.0
self.apply_params()
def gen_pulse(self, num_comp_waves=None,
min_wavelen=None, max_wavelen=None):
"""
Generate a random pulse by summing sine waves with random freq,
amplitude and phase offset
"""
if num_comp_waves is not None:
self.num_comp_waves = num_comp_waves
if min_wavelen is not None:
self.min_wavelen = min_wavelen
if max_wavelen is not None:
self.max_wavelen = max_wavelen
num_comp_waves = self.num_comp_waves
min_wavelen = self.min_wavelen
max_wavelen = self.max_wavelen
if min_wavelen > self.pulse_time:
raise ValueError("Minimum wavelength cannot be greater than "
"the pulse time")
if max_wavelen <= min_wavelen:
raise ValueError("Maximum wavelength must be greater than "
"the minimum wavelength")
if not self._pulse_initialised:
self.init_pulse()
# use some phase to avoid the first pulse being always 0
sum_wave = np.zeros(self.tau.shape)
t = self.time
wl_range = max_wavelen - min_wavelen
amp_scale = np.sqrt(8)*self.scaling / float(num_comp_waves)
for n in range(num_comp_waves):
amp = amp_scale*(np.random.rand()*2 - 1)
phase_off = np.random.rand()*np.pi/2.0
wavelen = min_wavelen + np.random.rand()*wl_range
curr_wave = amp*np.sin(2*np.pi*t/wavelen + phase_off)
sum_wave += curr_wave
return self._apply_bounds_and_offset(sum_wave)
class PulseGenRndWalk1(PulseGen):
"""
Generates pulses by using a random walk algorithm
Attributes
----------
scaling : float
Used as the range for the starting amplitude
Note must used bounds if values must be restricted.
Also scales the max_d_amp value
(copied from Dynamics.initial_ctrl_scaling if given)
max_d_amp : float
Maximum amount amplitude will change between timeslots
Note this is also factored by the scaling attribute
"""
def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.random = True
self.max_d_amp = 0.1
self.apply_params()
def gen_pulse(self, max_d_amp=None):
"""
Generate a pulse by changing the amplitude a random amount between
-max_d_amp and +max_d_amp at each timeslot. The walk will start at
a random amplitude between -/+scaling.
"""
if max_d_amp is not None:
self.max_d_amp = max_d_amp
max_d_amp = self.max_d_amp*self.scaling
if not self._pulse_initialised:
self.init_pulse()
walk = np.zeros(self.tau.shape)
amp = self.scaling*(np.random.rand()*2 - 1)
for k in range(len(walk)):
walk[k] = amp
amp += (np.random.rand()*2 - 1)*max_d_amp
return self._apply_bounds_and_offset(walk)
class PulseGenRndWalk2(PulseGen):
"""
Generates pulses by using a random walk algorithm
Note this is best used with bounds as the walks tend to wander far
Attributes
----------
scaling : float
Used as the range for the starting amplitude
Note must used bounds if values must be restricted.
Also scales the max_d2_amp value
(copied from Dynamics.initial_ctrl_scaling if given)
max_d2_amp : float
Maximum amount amplitude gradient will change between timeslots
Note this is also factored by the scaling attribute
"""
def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.random = True
self.max_d2_amp = 0.01
self.apply_params()
def gen_pulse(self, init_grad_range=None, max_d2_amp=None):
"""
Generate a pulse by changing the amplitude gradient a random amount
between -max_d2_amp and +max_d2_amp at each timeslot.
The walk will start at a random amplitude between -/+scaling.
The gradient will start at 0
"""
if max_d2_amp is not None:
self.max_d2_amp = max_d2_amp
max_d2_amp = self.max_d2_amp
if not self._pulse_initialised:
self.init_pulse()
walk = np.zeros(self.tau.shape)
amp = self.scaling*(np.random.rand()*2 - 1)
print("Start amp {}".format(amp))
grad = 0.0
print("Start grad {}".format(grad))
for k in range(len(walk)):
walk[k] = amp
grad += (np.random.rand()*2 - 1)*max_d2_amp
amp += grad
# print("grad {}".format(grad))
return self._apply_bounds_and_offset(walk)
[docs]class PulseGenLinear(PulseGen):
"""
Generates linear pulses
Attributes
----------
gradient : float
Gradient of the line.
Note this is calculated from the start_val and end_val if these
are given
start_val : float
Start point of the line. That is the starting amplitude
end_val : float
End point of the line.
That is the amplitude at the start of the last timeslot
"""
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.gradient = None
self.start_val = -1.0
self.end_val = 1.0
self.apply_params()
[docs] def init_pulse(self, gradient=None, start_val=None, end_val=None):
"""
Calculate the gradient if pulse is defined by start and
end point values
"""
PulseGen.init_pulse(self)
if start_val is not None and end_val is not None:
self.start_val = start_val
self.end_val = end_val
if self.start_val is not None and self.end_val is not None:
self.gradient = float(self.end_val - self.start_val) / \
(self.pulse_time - self.tau[-1])
[docs] def gen_pulse(self, gradient=None, start_val=None, end_val=None):
"""
Generate a linear pulse using either the gradient and start value
or using the end point to calulate the gradient
Note that the scaling and offset parameters are still applied,
so unless these values are the default 1.0 and 0.0, then the
actual gradient etc will be different
Returns the pulse as an array of vales for each timeslot
"""
if (gradient is not None or
start_val is not None or end_val is not None):
self.init_pulse(gradient, start_val, end_val)
if not self._pulse_initialised:
self.init_pulse()
pulse = np.empty(self.num_tslots)
t = 0.0
for k in range(self.num_tslots):
y = self.gradient*t + self.start_val
pulse[k] = self.scaling*y
t = t + self.tau[k]
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenPeriodic(PulseGen):
"""
Intermediate class for all periodic pulse generators
All of the periodic pulses range from -1 to 1
All have a start phase that can be set between 0 and 2pi
Attributes
----------
num_waves : float
Number of complete waves (cycles) that occur in the pulse.
wavelen and freq calculated from this if it is given
wavelen : float
Wavelength of the pulse (assuming the speed is 1)
freq is calculated from this if it is given
freq : float
Frequency of the pulse
start_phase : float
Phase of the pulse signal when t=0
"""
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.periodic = True
self.num_waves = None
self.freq = 1.0
self.wavelen = None
self.start_phase = 0.0
self.apply_params()
[docs] def init_pulse(self, num_waves=None, wavelen=None,
freq=None, start_phase=None):
"""
Calculate the wavelength, frequency, number of waves etc
from the each other and the other parameters
If num_waves is given then the other parameters are worked from this
Otherwise if the wavelength is given then it is the driver
Otherwise the frequency is used to calculate wavelength and num_waves
"""
PulseGen.init_pulse(self)
if start_phase is not None:
self.start_phase = start_phase
if num_waves is not None or wavelen is not None or freq is not None:
self.num_waves = num_waves
self.wavelen = wavelen
self.freq = freq
if self.num_waves is not None:
self.freq = float(self.num_waves) / self.pulse_time
self.wavelen = 1.0/self.freq
elif self.wavelen is not None:
self.freq = 1.0/self.wavelen
self.num_waves = self.wavelen*self.pulse_time
else:
self.wavelen = 1.0/self.freq
self.num_waves = self.wavelen*self.pulse_time
[docs]class PulseGenSine(PulseGenPeriodic):
"""
Generates sine wave pulses
"""
[docs] def gen_pulse(self, num_waves=None, wavelen=None,
freq=None, start_phase=None):
"""
Generate a sine wave pulse
If no params are provided then the class object attributes are used.
If they are provided, then these will reinitialise the object attribs.
returns the pulse as an array of vales for each timeslot
"""
if start_phase is not None:
self.start_phase = start_phase
if num_waves is not None or wavelen is not None or freq is not None:
self.init_pulse(num_waves, wavelen, freq, start_phase)
if not self._pulse_initialised:
self.init_pulse()
pulse = np.empty(self.num_tslots)
t = 0.0
for k in range(self.num_tslots):
phase = 2*np.pi*self.freq*t + self.start_phase
pulse[k] = self.scaling*np.sin(phase)
t = t + self.tau[k]
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenSquare(PulseGenPeriodic):
"""
Generates square wave pulses
"""
[docs] def gen_pulse(self, num_waves=None, wavelen=None,
freq=None, start_phase=None):
"""
Generate a square wave pulse
If no parameters are pavided then the class object attributes are used.
If they are provided, then these will reinitialise the object attribs
"""
if start_phase is not None:
self.start_phase = start_phase
if num_waves is not None or wavelen is not None or freq is not None:
self.init_pulse(num_waves, wavelen, freq, start_phase)
if not self._pulse_initialised:
self.init_pulse()
pulse = np.empty(self.num_tslots)
t = 0.0
for k in range(self.num_tslots):
phase = 2*np.pi*self.freq*t + self.start_phase
x = phase/(2*np.pi)
y = 4*np.floor(x) - 2*np.floor(2*x) + 1
pulse[k] = self.scaling*y
t = t + self.tau[k]
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenSaw(PulseGenPeriodic):
"""
Generates saw tooth wave pulses
"""
[docs] def gen_pulse(self, num_waves=None, wavelen=None,
freq=None, start_phase=None):
"""
Generate a saw tooth wave pulse
If no parameters are pavided then the class object attributes are used.
If they are provided, then these will reinitialise the object attribs
"""
if start_phase is not None:
self.start_phase = start_phase
if num_waves is not None or wavelen is not None or freq is not None:
self.init_pulse(num_waves, wavelen, freq, start_phase)
if not self._pulse_initialised:
self.init_pulse()
pulse = np.empty(self.num_tslots)
t = 0.0
for k in range(self.num_tslots):
phase = 2*np.pi*self.freq*t + self.start_phase
x = phase/(2*np.pi)
y = 2*(x - np.floor(0.5 + x))
pulse[k] = self.scaling*y
t = t + self.tau[k]
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenTriangle(PulseGenPeriodic):
"""
Generates triangular wave pulses
"""
[docs] def gen_pulse(self, num_waves=None, wavelen=None,
freq=None, start_phase=None):
"""
Generate a sine wave pulse
If no parameters are pavided then the class object attributes are used.
If they are provided, then these will reinitialise the object attribs
"""
if start_phase is not None:
self.start_phase = start_phase
if num_waves is not None or wavelen is not None or freq is not None:
self.init_pulse(num_waves, wavelen, freq, start_phase)
if not self._pulse_initialised:
self.init_pulse()
pulse = np.empty(self.num_tslots)
t = 0.0
for k in range(self.num_tslots):
phase = 2*np.pi*self.freq*t + self.start_phase + np.pi/2.0
x = phase/(2*np.pi)
y = 2*np.abs(2*(x - np.floor(0.5 + x))) - 1
pulse[k] = self.scaling*y
t = t + self.tau[k]
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenGaussian(PulseGen):
"""
Generates pulses with a Gaussian profile
"""
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self._uses_time = True
self.mean = 0.5*self.pulse_time
self.variance = 0.5*self.pulse_time
self.apply_params()
[docs] def gen_pulse(self, mean=None, variance=None):
"""
Generate a pulse with Gaussian shape. The peak is centre around the
mean and the variance determines the breadth
The scaling and offset attributes are applied as an amplitude
and fixed linear offset. Note that the maximum amplitude will be
scaling + offset.
"""
if not self._pulse_initialised:
self.init_pulse()
if mean:
Tm = mean
else:
Tm = self.mean
if variance:
Tv = variance
else:
Tv = self.variance
t = self.time
T = self.pulse_time
pulse = self.scaling*np.exp(-(t-Tm)**2/(2*Tv))
return self._apply_bounds_and_offset(pulse)
[docs]class PulseGenGaussianEdge(PulseGen):
"""
Generate pulses with inverted Gaussian ramping in and out
It's intended use for a ramping modulation, which is often required in
experimental setups.
Attributes
----------
decay_time : float
Determines the ramping rate. It is approximately the time
required to bring the pulse to full amplitude
It is set to 1/10 of the pulse time by default
"""
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self._uses_time = True
self.decay_time = self.pulse_time / 10.0
self.apply_params()
[docs] def gen_pulse(self, decay_time=None):
"""
Generate a pulse that starts and ends at zero and 1.0 in between
then apply scaling and offset
The tailing in and out is an inverted Gaussian shape
"""
if not self._pulse_initialised:
self.init_pulse()
t = self.time
if decay_time:
Td = decay_time
else:
Td = self.decay_time
T = self.pulse_time
pulse = 1.0 - np.exp(-t**2/Td) - np.exp(-(t-T)**2/Td)
pulse = pulse*self.scaling
return self._apply_bounds_and_offset(pulse)
### The following are pulse generators for the CRAB algorithm ###
# AJGP 2015-05-14:
# The intention is to have a more general base class that allows
# setting of general basis functions
[docs]class PulseGenCrab(PulseGen):
"""
Base class for all CRAB pulse generators
Note these are more involved in the optimisation process as they are
used to produce piecewise control amplitudes each time new optimisation
parameters are tried
Attributes
----------
num_coeffs : integer
Number of coefficients used for each basis function
num_basis_funcs : integer
Number of basis functions
In this case set at 2 and should not be changed
coeffs : float array[num_coeffs, num_basis_funcs]
The basis coefficient values
randomize_coeffs : bool
If True (default) then the coefficients are set to some random values
when initialised, otherwise they will all be equal to self.scaling
"""
def __init__(self, dyn=None, num_coeffs=None, params=None):
self.parent = dyn
self.num_coeffs = num_coeffs
self.params = params
self.reset()
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGen.reset(self)
self.NUM_COEFFS_WARN_LVL = 20
self.DEF_NUM_COEFFS = 4
self._BSC_ALL = 1
self._BSC_GT_MEAN = 2
self._BSC_LT_MEAN = 3
self._uses_time = True
self.time = None
self.num_basis_funcs = 2
self.num_optim_vars = 0
self.coeffs = None
self.randomize_coeffs = True
self._num_coeffs_estimated = False
self.guess_pulse_action = 'MODULATE'
self.guess_pulse = None
self.guess_pulse_func = None
self.apply_params()
[docs] def init_pulse(self, num_coeffs=None):
"""
Set the initial freq and coefficient values
"""
PulseGen.init_pulse(self)
self.init_coeffs(num_coeffs=num_coeffs)
if self.guess_pulse is not None:
self.init_guess_pulse()
self._init_bounds()
if self.log_level <= logging.DEBUG and not self._num_coeffs_estimated:
logger.debug(
"CRAB pulse initialised with {} coefficients per basis "
"function, which means a total of {} "
"optimisation variables for this pulse".format(
self.num_coeffs, self.num_optim_vars))
# def generate_guess_pulse(self)
# if isinstance(self.guess_pulsegen, PulseGen):
# self.guess_pulse = self.guess_pulsegen.gen_pulse()
# return self.guess_pulse
[docs] def init_coeffs(self, num_coeffs=None):
"""
Generate the initial ceofficent values.
Parameters
----------
num_coeffs : integer
Number of coefficients used for each basis function
If given this overides the default and sets the attribute
of the same name.
"""
if num_coeffs:
self.num_coeffs = num_coeffs
self._num_coeffs_estimated = False
if not self.num_coeffs:
if isinstance(self.parent, dynamics.Dynamics):
dim = self.parent.get_drift_dim()
self.num_coeffs = self.estimate_num_coeffs(dim)
self._num_coeffs_estimated = True
else:
self.num_coeffs = self.DEF_NUM_COEFFS
self.num_optim_vars = self.num_coeffs*self.num_basis_funcs
if self._num_coeffs_estimated:
if self.log_level <= logging.INFO:
logger.info(
"The number of CRAB coefficients per basis function "
"has been estimated as {}, which means a total of {} "
"optimisation variables for this pulse. Based on the "
"dimension ({}) of the system".format(
self.num_coeffs, self.num_optim_vars, dim))
# Issue warning if beyond the recommended level
if self.log_level <= logging.WARN:
if self.num_coeffs > self.NUM_COEFFS_WARN_LVL:
logger.warn(
"The estimated number of coefficients {} exceeds "
"the amount ({}) recommended for efficient "
"optimisation. You can set this level explicitly "
"to suppress this message.".format(
self.num_coeffs, self.NUM_COEFFS_WARN_LVL))
if self.randomize_coeffs:
r = np.random.random([self.num_coeffs, self.num_basis_funcs])
self.coeffs = (2*r - 1.0) * self.scaling
else:
self.coeffs = np.ones([self.num_coeffs,
self.num_basis_funcs])*self.scaling
[docs] def estimate_num_coeffs(self, dim):
"""
Estimate the number coefficients based on the dimensionality of the
system.
Returns
-------
num_coeffs : int
estimated number of coefficients
"""
num_coeffs = max(2, dim - 1)
return num_coeffs
[docs] def get_optim_var_vals(self):
"""
Get the parameter values to be optimised
Returns
-------
list (or 1d array) of floats
"""
return self.coeffs.ravel().tolist()
[docs] def set_optim_var_vals(self, param_vals):
"""
Set the values of the any of the pulse generation parameters
based on new values from the optimisation method
Typically this will be the basis coefficients
"""
# Type and size checking avoided here as this is in the
# main optmisation call sequence
self.set_coeffs(param_vals)
def set_coeffs(self, param_vals):
self.coeffs = param_vals.reshape(
[self.num_coeffs, self.num_basis_funcs])
def init_guess_pulse(self):
self.guess_pulse_func = None
if not self.guess_pulse_action:
logger.WARN("No guess pulse action given, hence ignored.")
elif self.guess_pulse_action.upper() == 'MODULATE':
self.guess_pulse_func = self.guess_pulse_modulate
elif self.guess_pulse_action.upper() == 'ADD':
self.guess_pulse_func = self.guess_pulse_add
else:
logger.WARN("No option for guess pulse action '{}' "
", hence ignored.".format(self.guess_pulse_action))
def guess_pulse_add(self, pulse):
pulse = pulse + self.guess_pulse
return pulse
def guess_pulse_modulate(self, pulse):
pulse = (1.0 + pulse)*self.guess_pulse
return pulse
def _init_bounds(self):
add_guess_pulse_scale = False
if self.lbound is None and self.ubound is None:
# no bounds to apply
self._bound_scale_cond = None
elif self.lbound is None:
# only upper bound
if self.ubound > 0:
self._bound_mean = 0.0
self._bound_scale = self.ubound
else:
add_guess_pulse_scale = True
self._bound_scale = self.scaling*self.num_coeffs + \
self.get_guess_pulse_scale()
self._bound_mean = -abs(self._bound_scale) + self.ubound
self._bound_scale_cond = self._BSC_GT_MEAN
elif self.ubound is None:
# only lower bound
if self.lbound < 0:
self._bound_mean = 0.0
self._bound_scale = abs(self.lbound)
else:
self._bound_scale = self.scaling*self.num_coeffs + \
self.get_guess_pulse_scale()
self._bound_mean = abs(self._bound_scale) + self.lbound
self._bound_scale_cond = self._BSC_LT_MEAN
else:
# lower and upper bounds
self._bound_mean = 0.5*(self.ubound + self.lbound)
self._bound_scale = 0.5*(self.ubound - self.lbound)
self._bound_scale_cond = self._BSC_ALL
def get_guess_pulse_scale(self):
scale = 0.0
if self.guess_pulse is not None:
scale = max(np.amax(self.guess_pulse) - np.amin(self.guess_pulse),
np.amax(self.guess_pulse))
return scale
def _apply_bounds(self, pulse):
"""
Scaling the amplitudes using the tanh function if there are bounds
"""
if self._bound_scale_cond == self._BSC_ALL:
pulse = np.tanh(pulse)*self._bound_scale + self._bound_mean
return pulse
elif self._bound_scale_cond == self._BSC_GT_MEAN:
scale_where = pulse > self._bound_mean
pulse[scale_where] = (np.tanh(pulse[scale_where])*self._bound_scale
+ self._bound_mean)
return pulse
elif self._bound_scale_cond == self._BSC_LT_MEAN:
scale_where = pulse < self._bound_mean
pulse[scale_where] = (np.tanh(pulse[scale_where])*self._bound_scale
+ self._bound_mean)
return pulse
else:
return pulse
[docs]class PulseGenCrabFourier(PulseGenCrab):
"""
Generates a pulse using the Fourier basis functions, i.e. sin and cos
Attributes
----------
freqs : float array[num_coeffs]
Frequencies for the basis functions
randomize_freqs : bool
If True (default) the some random offset is applied to the frequencies
"""
[docs] def reset(self):
"""
reset attributes to default values
"""
PulseGenCrab.reset(self)
self.freqs = None
self.randomize_freqs = True
[docs] def init_pulse(self, num_coeffs=None):
"""
Set the initial freq and coefficient values
"""
PulseGenCrab.init_pulse(self)
self.init_freqs()
[docs] def init_freqs(self):
"""
Generate the frequencies
These are the Fourier harmonics with a uniformly distributed
random offset
"""
self.freqs = np.empty(self.num_coeffs)
ff = 2*np.pi / self.pulse_time
for i in range(self.num_coeffs):
self.freqs[i] = ff*(i + 1)
if self.randomize_freqs:
self.freqs += np.random.random(self.num_coeffs) - 0.5
[docs] def gen_pulse(self, coeffs=None):
"""
Generate a pulse using the Fourier basis with the freqs and
coeffs attributes.
Parameters
----------
coeffs : float array[num_coeffs, num_basis_funcs]
The basis coefficient values
If given this overides the default and sets the attribute
of the same name.
"""
if coeffs:
self.coeffs = coeffs
if not self._pulse_initialised:
self.init_pulse()
pulse = np.zeros(self.num_tslots)
for i in range(self.num_coeffs):
phase = self.freqs[i]*self.time
# basis1comp = self.coeffs[i, 0]*np.sin(phase)
# basis2comp = self.coeffs[i, 1]*np.cos(phase)
# pulse += basis1comp + basis2comp
pulse += self.coeffs[i, 0]*np.sin(phase) + \
self.coeffs[i, 1]*np.cos(phase)
if self.guess_pulse_func:
pulse = self.guess_pulse_func(pulse)
if self.ramping_pulse is not None:
pulse = self._apply_ramping_pulse(pulse)
return self._apply_bounds(pulse)