Source code for qubiter.StateVec

import pprint as pp
import scipy as sc
import qubiter.utilities_gen as ut
from collections import OrderedDict
import pandas as pan
from qubiter.Plotter import *

import sys
if 'autograd.numpy' not in sys.modules:
    import numpy as np
    import autograd.numpy as np

[docs]class StateVec: """ This class is a wrapper for its main attribute, a complex numpy array self.arr with shape [2]*num_qbits. The class also provides functions for performing calculations on dictionaries of the objects of this class StateVec. The keys of these dictionaries of state vectors are strings that we call branch_keys, because they name "branches" in class SEO_simulation. This class also provides a function for constructing from such dictionaries of state vectors, a density matrix which is a 2 dim square numpy array of dimension 2^num_qbits. IMPORTANT: See docstring of method get_traditional_st_vec() for explanation of qubit ordering conventions and shape of self.arr. Attributes ---------- arr : np.ndarray a complex array of shape [2]*num_qbits num_qbits : int """
[docs] def __init__(self, num_qbits, arr=None): """ Constructor Parameters ---------- num_qbits : int arr : np.ndarray Returns ------- """ self.num_qbits = num_qbits self.arr = arr if arr is not None: assert self.arr.shape == tuple([2]*self.num_qbits)
[docs] @staticmethod def is_zero(st_vec): """ Returns True iff an object of this class is None or its parameter `arr` is None. Parameters ---------- st_vec : StateVec|None Returns ------- bool """ return st_vec is None or st_vec.arr is None
[docs] def __str__(self): """ Returns str(self.arr). Returns ------- str """ return str(self.arr)
[docs] @staticmethod def get_ground_st_vec(num_qbits): """ Returns StateVec for the ground state ``|0>|0>|0>...|0>``, where ``|0> = [ 1,0]^t`` and ``|1> = [0,1]^t``, t = transpose. Parameters ---------- num_qbits : int Returns ------- StateVec """ ty = np.complex128 assert num_qbits > 0 arr = np.zeros([1 << num_qbits], dtype=ty) arr[0] = 1 arr = arr.reshape([2]*num_qbits) return StateVec(num_qbits, arr)
[docs] @staticmethod def get_random_st_vec(num_qbits, rand_seed=None): """ Returns StateVec for random state ``\sum_b^n A(b^n)|b^n>``, ``b^n \in { 0, 1}^n``, where n=num_qbits and ``\sum_b^n |A( b^n)|^2 = 1``. Parameters ---------- num_qbits : int rand_seed : int Returns ------- StateVec """ if rand_seed: np.random.seed(rand_seed) # returns array of random numbers in [0, 1] interval mat_phi = 2*np.pi*np.random.random(1 << num_qbits) mat_r = np.random.random(1 << num_qbits) arr = mat_r*(np.cos(mat_phi) + 1j*np.sin(mat_phi)) magnitude = np.linalg.norm(arr) arr /= magnitude arr = arr.reshape([2]*num_qbits) return StateVec(num_qbits, arr)
[docs] @staticmethod def get_standard_basis_st_vec(spin_dir_list, ZL=True): """ If ZL = True, returns StateVec for state ``...|s2>|s1>|s0>``, where ``spin_dir_list=[...,s2, s1, s0], s_j \in {0, 1}`` for all j, ``|0> = [1, 0]^t`` and ``|1> = [0,1]^t``, t = transpose. If ZL = False, same except ``spin_dir_list=reversed([...,s2, s1, s0])``. Parameters ---------- spin_dir_list : list[int] ZL : bool True(False) if last(first) entry of spin_dir_list refers to qubit 0 Returns ------- StateVec """ num_qbits = len(spin_dir_list) arr = np.zeros([1 << num_qbits], dtype=np.complex128) arr = arr.reshape([2]*num_qbits) if ZL: spin_dir_list = reversed(spin_dir_list) # print("spins", list(spin_dir_list)) # print("arr", arr.shape) arr[tuple(spin_dir_list)] = 1 return StateVec(num_qbits, arr)
@property def get_traditional_st_vec(self): """ **IMPORTANT: Internally, self.arr in Qubiter has shape [2]*num_qbits and assumes ZF convention because that way a numpy axis and a qubit number are the same thing. However, the traditional way of writing a state vector is as a column array of dimension 1<< num_qbits in the ZL convention.** This function returns the traditional view. So it reshapes ( flattens) the array and it reverses the axes (reversing axes takes it from ZF to ZL). The rows are always labelled 0, 1, 2, 3, ... or the binary representation thereof, regardless of whether ZL or ZF convention. One can go from digital to binary labels and vice versa using >>> x = np.binary_repr(3, width=4) >>> x '0011' >>> int(x, 2) 3 Parameters ---------- Returns ------- np.array """ perm = list(reversed(range(self.num_qbits))) return np.transpose(self.arr, perm).flatten()
[docs] @staticmethod def get_den_mat(num_qbits, st_vec_dict): """ Returns a density matrix (indexed in ZL convention) constructed from st_vec_dict which is a dict from strings to StateVec. The rows and columns are always labelled 0, 1, 2, .. or binary representation thereof, regardless of whether ZL or ZF convention. To switch between bin to dec representations of labels, see docstring of get_traditional_st_vec(). Parameters ---------- num_qbits : int st_vec_dict : dict[str, StateVec] Returns ------- np.ndarray """ dim = 1 << num_qbits den_mat = np.zeros((dim, dim), dtype=complex) # print(",,,", den_mat) for br_key in st_vec_dict: if StateVec.is_zero(st_vec_dict[br_key]): continue vec = st_vec_dict[br_key].get_traditional_st_vec assert vec.shape == (dim,) den_mat += np.outer(vec, np.conj(vec)) # print(',,..', br_key, vec) tr = np.trace(den_mat) assert abs(tr) > 1e-6 return den_mat/tr
[docs] @staticmethod def get_partial_tr(num_qbits, den_mat, traced_bits_set): """ Returns the partial trace of a density matrix den_mat. Traces over qubits in set traced_bits_set. To get full trace, just do np.trace( den_mat). Parameters ---------- num_qbits : int den_mat : np.ndarray if dim=2^num_qbits, this function assumes that den_mat has shape (dim, dim) and that it's indexed in the ZL convention so qubit 0 corresponds to axis num_qbits-1. traced_bits_set : set[int] Set of qubits being traced over Returns ------- np.ndarray """ dim = 1 << num_qbits assert den_mat.shape == (dim, dim) assert set(range(num_qbits)) > traced_bits_set dm = den_mat.reshape([2]*(2*num_qbits)) # bit 0 corresponds to axis num_qbits - 1 traced_axes = [num_qbits - 1 - k for k in traced_bits_set] num_traces = len(traced_axes) for k in range(num_traces): ax = traced_axes.pop(0) dm = np.trace(dm, axis1=ax, axis2=ax + num_qbits - k) traced_axes = list(map(lambda x: (x if x <= ax else x-1), traced_axes)) new_num_qbits = num_qbits - len(traced_bits_set) dim = 1 << new_num_qbits dm = dm.reshape((dim, dim)) return dm
[docs] @staticmethod def get_impurity(den_mat): """ Returns abs(trace(den_mat^2) -1). This is zero iff the density matrix den_mat represents a pure state. For example, for a pure state ``den_mat = |a><a|``, ``den_mat^2 = den_mat = |a><a|`` so this quantity is indeed zero. Parameters ---------- den_mat : np.nparray density matrix, shape=(dim, dim) where dim=2^num_qbits Returns ------- float """ return abs(np.trace(, den_mat)) - 1)
[docs] @staticmethod def get_entropy(den_mat, method='eigen'): """ Returns entropy of density matrix den_mat. Uses natural log for entropy. Parameters ---------- den_mat : np.ndarray Density matrix. Eigenvalues must be non-negative and sum to 1. method : str method used to calculate log of array. Either 'eigen' or 'pade'. Returns ------- float """ ent = 0.0 if method == 'eigen': evas = np.real(np.linalg.eigvalsh(den_mat)) assert np.all(evas > -1e-6), evas assert abs(np.sum(evas) - 1) < 1e-6 for val in evas: if val > 1e-6: ent += - val*np.log(val) elif method == 'pade': ent = - np.trace(, sc.linalg.logm(den_mat))) else: assert False, 'unsupported method for ' + \ 'calculating entropy of a density matrix.' return ent
[docs] @staticmethod def get_den_mat_pd(den_mat): """ Returns the diagonal of den_mat (so indexed in ZL convention). den_mat is expected to be a density matrix returned by get_den_mat(). Parameters ---------- den_mat : np.ndarray density matrix, shape=(dim, dim) where dim=2^num_qbits, indexed in ZL convention. Returns ------- np.ndarray """ return np.real(np.diag(den_mat))
[docs] def get_pd(self): """ Returns copy of self.get_traditional_st_vec() with amplitudes replaced by probabilities. pd = probability distribution. So returns one column array indexed in ZL convention like the traditional state vec is. Doesn't check that the resulting array sums to 1. Parameters ---------- Returns ------- np.ndarray probability distribution of shape (2^num_qbits,) IMP: will be indexed in ZL convention. """ x = self.get_traditional_st_vec return np.real(x*np.conj(x))
[docs] def get_mean_value_of_real_diag_mat(self, real_arr): """ In Quantum Mechanics, one often needs to calculate the mean value of a Hermitian operator H, ``mean= <psi|H|psi>``. Let ``H = U^\dag D U``, where U is unitary and D is real diagonal matrix. If ``self = U|psi>``, then this reduces to finding ``mean= <self|D|self>``. So must decompose U into a SEO and evolve, using SEO_simulator, to the state ``U|psi>``. Parameters ---------- real_arr : np.ndarray a real array of shape=[2]^num_qbits (same shape as self.arr). If flattened, real_arr contains the diagonal of the matrix D. If U is a Kronecker prod of 2-dim unitary matrices, the flattened real_arr can be obtained as Kronecker product of spinors, i.e., shape=( 2, ) arrays. Returns ------- float """ return np.sum(np.real(np.conj(self.arr) * real_arr * self.arr))
[docs] def get_total_prob(self): """ Returns total probability of self. Parameters ---------- Returns ------- float """ return np.sum(np.real(np.conj(self.arr)*self.arr))
[docs] @staticmethod def get_observations_vec(num_qbits, pd, num_shots, rand_seed=None): """ vec = vector num_shots (number of shots) is often called number of trials or number of samples. For num_shots=1, this method returns an int (actually, a 1 X 1 array with an int in it) in `range(1<<num_qbits)` chosen according to the probability distribution pd for num_qbits qubits. If the output int were to be expressed in binary notation, its last, rightmost bit would be the measurement of the 0th qubit (because pd is assumed to be in ZL convention). For num_shots>1, the method returns an np.ndarray of shape ( num_shots,) with the result of doing num_shots repetitions of what was done for num_shots=1. Does not assume that pd is normalized to 1. Parameters ---------- num_qbits : int pd : np.ndarray probability distribution of shape (2^num_qbits,) IMP: assumed to be indexed in ZL convention num_shots : int rand_seed : int Returns ------- np.ndarray shape (num_shots,) """ if rand_seed: np.random.seed(rand_seed) len_pd = 1 << num_qbits assert pd.shape == (len_pd,) tot_prob = np.sum(pd) p = pd if abs(tot_prob-1) > 1e-5: p = pd/tot_prob return np.random.choice(np.arange(0, len_pd), size=num_shots, p=p)
[docs] @staticmethod def get_counts_from_obs_vec(num_qbits, obs_vec, use_bin_labels=True, omit_zero_counts=True): """ This method takes as input an observations vector obs_vec such as returned by another method in this class, namely get_observations_vec(). This method returns an OrderedDict called state_name_to_count that maps the names of states to the number of times they occur in obs_vec. If use_bin_labels=True, state names are a string composed of a binary number that is num_qbits long, followed by 'ZL' because ZL convention is assumed. If use_bin_labels=False, state names are '0', '1', '2', etc. Parameters ---------- num_qbits : int obs_vec : np.ndarray use_bin_labels : bool omit_zero_counts : bool Returns ------- OrderedDict[str, int] """ obs_list = list(obs_vec) num_states = 1 << num_qbits state_name_to_count = OrderedDict() for s in range(num_states): s_count = obs_list.count(s) if use_bin_labels: # this returns string key = np.binary_repr(s, width=num_qbits) key += 'ZL' else: key = str(s) if s_count > 0 or not omit_zero_counts: state_name_to_count[key] = s_count return state_name_to_count
[docs] @staticmethod def get_empirical_pd_from_counts(num_qbits, state_name_to_count): """ This method takes as input "the counts dict" (i.e., an OrderedDict called state_name_to_count which is produced by another method in this class, namely get_counts_from_obs_vec()). This method returns an empirical probability distribution emp_pd calculated from the counts dict. emp_pd indices are ints referring to qubit states labelled in the ZL convention. Parameters ---------- num_qbits : int state_name_to_count : OrderedDict[str, int] Returns ------- emp_pd : np.ndarray its shape is (1<<num_qbits,) """ emp_pd = np.zeros(shape=(1 << num_qbits,), dtype=float) tot_counts = 0 for st_name, count in state_name_to_count.items(): # state name ends in ZL so trim last two chars pos = int(st_name[:-2], 2) emp_pd[pos] = count tot_counts += count return emp_pd/tot_counts
[docs] @staticmethod def get_emp_state_vec_from_emp_pd(num_qbits, emp_pd): """ This method takes as input an empirical probability distribution emp_pd and it returns an empirical state vector calculated from emp_pd. This requires reshaping emp_pd to the shape [2]*num_qbits, permuting its indices from the ZL to the ZF convention, and then taking the sqrt of the components to get an amplitude instead of a probability. All amplitudes of the output state vector are real though. Parameters ---------- num_qbits : int emp_pd : np.ndarray its shape is (1<<num_qbits,) Returns ------- StateVec """ assert emp_pd.shape == (1 << num_qbits,) arr = emp_pd.reshape(tuple([2]*num_qbits)) perm = list(reversed(range(num_qbits))) arr = np.transpose(arr, perm) sqrt_probs = np.sqrt(arr) return StateVec(num_qbits, sqrt_probs)
[docs] @staticmethod def get_bit_probs(num_qbits, pd): """ Returns a list whose jth item is, for the jth qubit, the pair (p, 1-p), where p is the probability that the jth qubit is 0, if the state of all other qubits is ignored. Does not assume that pd is normalized to 1. Parameters ---------- num_qbits : int pd : np.ndarray probability distribution of shape (2^num_qbits,) IMP: assumed to be indexed in ZL convention Returns ------- list[tuple[float, float]] """ assert pd.shape == (1 << num_qbits,) probs = [] arr = pd.reshape([2] * num_qbits) # tot_prob may not be one # if a measurement has been done tot_prob = np.sum(arr) if abs(tot_prob-1) > 1e-5: arr /= tot_prob # slicex is a portmanteau of slice index # print("state_vec_pd=\n", vec) slicex = [slice(None)]*num_qbits # pd is assumed to be in ZL convention # so when reshape, bit k has axis= num_qbits -1 - k for k in range(num_qbits): k_axis = num_qbits - 1 - k slicex[k_axis] = 0 p = np.sum(arr[tuple(slicex)]) probs.append((p, 1-p)) slicex[k_axis] = slice(None) # restore to all entries slice(None) # print(probs) return probs
# this method is correct but superfluous, I think # @staticmethod # def get_bit_counts(bit_probs, num_shots, rand_seed=None): # """ # num_shots (number of shots) is often called number of trials or # number of samples. # # Returns a list whose jth item is, for the jth qubit, the pair ( # count0, count1), where count0 + count1 = num_shots, and as num_shots # ->infinity, (count0, count1)/num_shots tends to the probability pair # bit_probs[j]. # # Parameters # ---------- # bit_probs : list[tuple[float, float]] # num_shots : int # rand_seed : int # # Returns # ------- # list[tuple[int, int]] # # """ # if rand_seed: # np.random.seed(rand_seed) # num_qbits = len(bit_probs) # counts = [] # for bit in range(num_qbits): # x1 = np.random.binomial(n=num_shots, p=bit_probs[bit][1]) # x0 = num_shots - x1 # counts.append((x0, x1)) # return counts
[docs] def pp_arr_entries(self, omit_zero_amps=False, show_pp_probs=False, ZL=True): """ pp=pretty print. Prints for each entry of self.arr, a line of the form (i, j, k, ...) self.arr[i, j, k, ...], with zero bit last ( resp., first) if ZL=True (resp., False). Parameters ---------- omit_zero_amps : bool If True, will not list states with zero amplitude. show_pp_probs : bool If True, will show probability of each amplitude. ZL : bool If True, multi-index in ZL (Zero bit Last) convention. If False, multi-index in ZF (Zero bit First) convention. Returns ------- None """ for ind, x in np.ndenumerate(self.arr): index = ind label = 'ZF' if ZL: index = ind[::-1] # this reverses order of tuple label = 'ZL' # turn index tuple into string and remove commas ind_str = str(index).replace(', ', '') + label mag = np.absolute(x) extra_str = '' if show_pp_probs: extra_str = '\t prob=' + '{0:6f}'.format(mag**2) x_str = ' (' if x.real < 0 else ' ( ' x_str += '{0:.6f} {1} {2:.6f}j)'.\ format(x.real, '-' if x.imag < 0 else '+', abs(x.imag)) if omit_zero_amps: if mag > 1E-6: print(ind_str + x_str + extra_str) else: print(ind_str + x_str + extra_str)
[docs] @staticmethod def get_style_dict(style): """ Given a style string as input, this method returns a dict mapping various strings denoting parameters of the method StateVec.describe_self() to their bool values for the input style. Parameters ---------- style : str Returns ------- dict[str, bool] """ vanilla = { 'print_st_vec': False, 'do_pp': False, 'omit_zero_amps': False, 'show_pp_probs': False, 'ZL': True, 'plot_st_vec_pd': False} if style == 'V1': out = vanilla elif style == 'ALL': out = {x: True for x in vanilla.keys()} out['plot_st_vec_pd'] = False elif style == 'ALL+': out = {x: True for x in vanilla.keys()} else: assert False, "unsupported style for StateVec.describe_self()" return out
[docs] def describe_self(self, print_st_vec=False, do_pp=False, omit_zero_amps=False, show_pp_probs=False, ZL=True, plot_st_vec_pd=False): """ Prints a description of self. Parameters ---------- print_st_vec : bool if True, prints the final state vector (which may be huge. For n qubits, it has 2^n components.) do_pp : bool pp= pretty print. Only used if print_st_vec=True. For pp=False, it prints final state vector in usual numpy array print style. For pp=True, it prints final state vector as column of (index, array value) pairs. omit_zero_amps : bool If print_st_vec=True, pp=True and this parameter is True too, will omit states with zero amplitude show_pp_probs : bool If True, will show probability of each standard basis state. ZL : bool If True, multi-index of ket in ZL (Zero bit Last) convention. If False, multi-index of ket in ZF (Zero bit First) convention. plot_st_vec_pd : bool If True, plots state vector's probability distribution. Returns ------- None """ if self.arr is None: print("zero state vector") return if print_st_vec: print('state vector:') if do_pp: if not ZL: print('ZF convention (Zero bit First in state tuple)') else: print('ZL convention (Zero bit Last in state tuple)') self.pp_arr_entries(omit_zero_amps, show_pp_probs, ZL) else: print(self.arr) print('total probability of state vector ' + '(=one if no measurements)=', "{0:0.6f}".format(self.get_total_prob())) print('dictionary with key=qubit, value=(Prob(0), Prob(1))') bit_probs = StateVec.get_bit_probs(self.num_qbits, self.get_pd()) bit_probs = [(round(x, 6), round(y, 6)) for x, y in bit_probs] pp.pprint(dict(enumerate(bit_probs))) if plot_st_vec_pd: st_vec_pd = self.get_pd() st_vec_pd_df = Plotter.get_pd_df(self.num_qbits, st_vec_pd) Plotter.plot_probs_col(['st_vec_pd'], [st_vec_pd_df])
[docs] @staticmethod def describe_st_vec_dict(st_vec_dict, **kwargs): """ Calls describe_self() for each branch of st_vec_dict. Parameters ---------- st_vec_dict : dict[str, StateVec] kwargs : dict[] the keyword arguments of describe_self() Returns ------- None """ for br_key in st_vec_dict.keys(): print("*********branch= " + br_key) if StateVec.is_zero(st_vec_dict[br_key]): print("zero state vector") else: st_vec_dict[br_key].describe_self(**kwargs)
if __name__ == "__main__": def main(): num_qbits = 3 gs = StateVec(num_qbits, arr=StateVec.get_ground_st_vec(num_qbits).arr) print('gs=\n', gs.arr) print("gs_trad=\n", gs.get_traditional_st_vec) S0100_ZL = StateVec(4, arr=StateVec.get_standard_basis_st_vec([0, 1, 0, 0], ZL=True).arr) print("S0100_ZL=\n", S0100_ZL.get_traditional_st_vec) S0100_ZF = StateVec(4, arr=StateVec.get_standard_basis_st_vec([0, 1, 0, 0], ZL=False).arr) print("S0100_ZF=\n", S0100_ZF.get_traditional_st_vec) st_vec0 = StateVec(num_qbits, arr=StateVec.get_random_st_vec(num_qbits).arr) st_vec1 = StateVec(num_qbits, arr=StateVec.get_random_st_vec(num_qbits).arr) st_vec_dict = {'br0': st_vec0, 'br1': st_vec1, 'br3': None} trad_st_vec = st_vec0.get_traditional_st_vec den_mat = StateVec.get_den_mat(num_qbits, st_vec_dict) print("den_mat\n", den_mat) print('trace_02 den_mat\n', StateVec.get_partial_tr(num_qbits, den_mat, {0, 2})) print("impurity=", StateVec.get_impurity(den_mat)) print("entropy=", StateVec.get_entropy(den_mat)) den_mat_pd = StateVec.get_den_mat_pd(den_mat) print('den_mat_pd=', den_mat_pd) st_vec_pd = st_vec0.get_pd() bit_probs_vec = StateVec.get_bit_probs(num_qbits, st_vec_pd) bit_probs_dm = StateVec.get_bit_probs(num_qbits, den_mat_pd) # print("counts_dm=\n", StateVec.get_bit_counts(bit_probs_dm, 10)) obs_vec = StateVec.get_observations_vec( num_qbits, st_vec_pd, num_shots=20) print('observations vec\n', obs_vec) counts = StateVec.get_counts_from_obs_vec(num_qbits, obs_vec) print('counts from obs_vec\n', counts) StateVec.describe_st_vec_dict(st_vec_dict, print_st_vec=True, do_pp=True, omit_zero_amps=False, show_pp_probs=True, ZL=True) main() import doctest doctest.testmod(verbose=True)