Source code for qutip.control.pulsegen

# -*- 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.control.dynamics as dynamics
import qutip.control.errors as errors


[docs]def create_pulse_gen(pulse_type='RND', dyn=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) if pulse_type == 'RNDFOURIER': return PulseGenRndFourier(dyn) if pulse_type == 'RNDWAVES': return PulseGenRndWaves(dyn) if pulse_type == 'RNDWALK1': return PulseGenRndWalk1(dyn) if pulse_type == 'RNDWALK2': return PulseGenRndWalk2(dyn) elif pulse_type == 'LIN': return PulseGenLinear(dyn) elif pulse_type == 'ZERO': return PulseGenZero(dyn) elif pulse_type == 'SINE': return PulseGenSine(dyn) elif pulse_type == 'SQUARE': return PulseGenSquare(dyn) elif pulse_type == 'SAW': return PulseGenSaw(dyn) elif pulse_type == 'TRIANGLE': return PulseGenTriangle(dyn) else: raise ValueError("No option for pulse_type '{}'".format(pulse_type))
[docs]class PulseGen: """ 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 """ def __init__(self, dyn=None): self.parent = dyn 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 else: self.num_tslots = 100 self.pulse_time = 1.0 self.scaling = 1.0 self.tau = None self.offset = 0.0 self._pulse_initialised = False self.periodic = False self.random = False self.lbound = -np.Inf self.ubound = np.Inf
[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 self._pulse_initialised = True 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 np.isinf(self.lbound) and np.isinf(self.ubound): return pulse + self.offset max_amp = max(pulse) min_amp = min(pulse) if (max_amp + self.offset <= self.ubound and min_amp + self.offset >= self.lbound): return pulse + self.offset # Some shifting / scaling is required. bound_range = self.ubound - self.lbound if np.isinf(bound_range): # One of the bounds is inf, so just shift the pulse if np.isinf(self.lbound): # 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: 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)
[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
[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 try: self.min_wavelen = self.pulse_time / 10.0 except: self.min_wavelen = 0.1 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 = np.zeros(self.num_tslots, dtype=float) for k in range(self.num_tslots-1): t[k+1] = t[k] + self.tau[k] 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.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 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 = np.zeros(self.num_tslots, dtype=float) for k in range(self.num_tslots-1): t[k+1] = t[k] + self.tau[k] 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 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 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
[docs] def init_pulse(self, gradient=None, start_val=None, end_val=None): """ Calulate 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
[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 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)