The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def apply_gate(self,gate,on_qubit_name):
on_qubit=self.qubits.get_quantum_register_containing(on_qubit_name)
if len(on_qubit.get_noop()) > 0:
print "NOTE this qubit has been measured previously, there should be no more gates allowed but we are reverting that measurement for consistency with IBM's language"
on_qubit.set_state(on_qubit.get_noop())
on_qubit.set_noop([])
if not on_qubit.is_entangled():
if on_qubit.get_num_qubits()!=1:
raise Exception("This qubit is not marked as entangled but it has an entangled state")
on_qubit.set_state(gate*on_qubit.get_state())
else:
if not on_qubit.get_num_qubits()>1:
raise Exception("This qubit is marked as entangled but it does not have an entangled state")
n_entangled=len(on_qubit.get_entangled())
apply_gate_to_qubit_idx=[qb.name for qb in on_qubit.get_entangled()].index(on_qubit_name)
if apply_gate_to_qubit_idx==0:
entangled_gate=gate
else:
entangled_gate=Gate.eye
for i in range(1,n_entangled):
if apply_gate_to_qubit_idx==i:
entangled_gate=np.kron(entangled_gate,gate)
else:
entangled_gate=np.kron(entangled_gate,Gate.eye)
on_qubit.set_state(entangled_gate*on_qubit.get_state()) Example 2
def defineSensorLoc(self,XYZ):
#############################################
# DEFINE TRANSMITTER AND RECEIVER LOCATIONS
#
# XYZ: N X 3 array containing transmitter center locations
# **NOTE** for this sensor, we know where the receivers are relative to transmitters
self.TxLoc = XYZ
dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
dx = mkvc(dx)
dy = mkvc(dy)
N = np.shape(XYZ)[0]
X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
Z = np.kron(XYZ[:,2],np.ones((25)))
self.RxLoc = np.c_[X,Y,Z]
self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25))) Example 3
def defineSensorLoc(self,XYZ):
#############################################
# DEFINE TRANSMITTER AND RECEIVER LOCATIONS
#
# XYZ: N X 3 array containing transmitter center locations
# **NOTE** for this sensor, we know where the receivers are relative to transmitters
self.TxLoc = XYZ
dx = np.kron([-0.18,0.,0.,0.,0.18],np.ones(3))
dy = np.kron([0.,-0.18,0.,0.18,0.],np.ones(3))
N = np.shape(XYZ)[0]
X = np.kron(XYZ[:,0],np.ones((15))) + np.kron(np.ones((N)),dx)
Y = np.kron(XYZ[:,1],np.ones((15))) + np.kron(np.ones((N)),dy)
Z = np.kron(XYZ[:,2],np.ones((15)))
self.RxLoc = np.c_[X,Y,Z]
self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((15)))
self.RxComp = np.kron(np.kron(np.ones(np.shape(XYZ)[0]),np.ones((5))),np.r_[1,2,3]) Example 4
def toa_calc_d_from_xy(x, y):
dimx = x.shape
x_dim = dimx[0]
m = dimx[1]
dimy = y.shape
y_dim = dimy[0]
n = dimy[1]
if x_dim > y_dim:
y[(y_dim + 1):x_dim, ] = np.zeros(x_dim - y_dim, n)
elif y_dim > x_dim:
x[(x_dim + 1):y_dim, ] = np.zeros(y_dim - x_dim, n)
d = sqrt(
((np.kron(np.ones(1, n), x) - np.kron(y, np.ones(1, m))) ** 2).sum(
axis=0))
d = np.asarray(d)
d = np.reshape(d, (m, n), order='F')
return d Example 5
def edgeCurl(self):
"""The edgeCurl property."""
if self.nCy > 1:
raise NotImplementedError('Edge curl not yet implemented for '
'nCy > 1')
if getattr(self, '_edgeCurl', None) is None:
# 1D Difference matricies
dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1, 0],
self.nCx, self.nCx, format="csr")
dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0, 1],
self.nCz, self.nCz+1, format="csr")
# 2D Difference matricies
Dr = sp.kron(sp.identity(self.nNz), dr)
Dz = -sp.kron(dz, sp.identity(self.nCx))
A = self.area
E = self.edge
# Edge curl operator
self._edgeCurl = (utils.sdiag(1/A)*sp.vstack((Dz, Dr)) *
utils.sdiag(E))
return self._edgeCurl Example 6
def aveE2CC(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CC', None) is None:
# The number of cell centers in each direction
n = self.vnC
if self.isSymmetric:
avR = utils.av(n[0])[:, 1:]
avR[0, 0] = 1.
self._aveE2CC = sp.kron(utils.av(n[2]), avR, format="csr")
else:
raise NotImplementedError('wrapping in the averaging is not '
'yet implemented')
# self._aveE2CC = (1./3)*sp.hstack((utils.kron3(utils.av(n[2]),
# utils.av(n[1]),
# utils.speye(n[0])),
# utils.kron3(utils.av(n[2]),
# utils.speye(n[1]),
# utils.av(n[0])),
# utils.kron3(utils.speye(n[2]),
# utils.av(n[1]),
# utils.av(n[0]))),
# format="csr")
return self._aveE2CC Example 7
def f(W,G,y,hparams):
# f = -1/N*sum_t log(exp(w(yt)'gt)/sum_k exp(wk'gt)) + l*||W||
# = -1/N*sum_t [w(yt)'*gt - log(sum_k exp(wk'gt))] + l*||W||
# = -1/N*sum(sum(W(:,y).*G,1),2) + 1/N*sum(log(sumexpWG),2) + l*sum(sum(W.^2));
#K,l = hparams
K = hparams['K']
l = hparams['l']
d,N = G.shape
W = W.reshape((d,K))
WG = np.dot(W.T,G) # K x N
WG -= np.kron(np.ones((K,1)),WG.max(axis=0).reshape(1,N))
#WG_max = WG.max(axis=0).reshape((1,N))
#expWG = np.exp(WG-np.kron(np.ones((K,1)),WG_max)) # K x N
expWG = np.exp(WG) # K x N
sumexpWG = expWG.sum(axis=0) # N x 1
WyG = WG[y,range(N)]
#WyG -= WG_max
fval = -1.0/N*(WyG).sum() \
+ 1.0/N*np.log(sumexpWG).sum() \
+ l*(W**2).sum()#(axis=(0,1))
return fval Example 8
def dfdv(W,G,y,hparams):
# df/dwk = -1/N*sum(x(:,y==k),2) + 1/N*sum_t exp(wk'xt)*xt/(sum_k exp(wk'xt))] + l*2*wk
K = hparams['K']
l = hparams['l']
d,N = G.shape
shapeW = W.shape
W = W.reshape((d,K))
WG = np.dot(W.T,G) # K x N
WG -= np.kron(np.ones((K,1)),WG.max(axis=0).reshape(1,N))
expWG = np.exp(WG) # K x N
sumexpWG = expWG.sum(axis=0) # N x 1
df = np.zeros((d,K))
for k in range(K):
indk = np.where(y==k)[0]
df[:,k] = -1./N*G[:,indk].sum(axis=1).reshape((d,)) \
+ 1./N*np.dot(G,(expWG[k,:]/sumexpWG).T).reshape((d,)) \
+ 2.*l*W[:,k].reshape((d,))
assert np.isnan(df).any()==False
return df.reshape(shapeW) Example 9
def estimator_cov(self,method):
""" Creates covariance matrix for the estimators
Parameters
----------
method : str
Estimation method
Returns
----------
A Covariance Matrix
"""
Y = np.array([reg[self.lags:] for reg in self.data])
Z = self._create_Z(Y)
if method == 'OLS':
sigma = self.ols_covariance()
else:
sigma = self.custom_covariance(self.latent_variables.get_z_values())
return np.kron(np.linalg.inv(np.dot(Z,np.transpose(Z))), sigma) Example 10
def entangling_mat(gate):
"""
Helper function to create the entangling gate matrix
"""
echoCR = expm(1j * pi / 4 * np.kron(pX, pZ))
if gate == "CNOT":
return echoCR
elif gate == "iSWAP":
return reduce(lambda x, y: np.dot(y, x),
[echoCR, np.kron(C1[6], C1[6]), echoCR])
elif gate == "SWAP":
return reduce(lambda x, y: np.dot(y, x),
[echoCR, np.kron(C1[6], C1[6]), echoCR, np.kron(
np.dot(C1[6], C1[1]), C1[1]), echoCR])
else:
raise ValueError("Entangling gate must be one of: CNOT, iSWAP, SWAP.") Example 11
def clifford_mat(c, numQubits):
"""
Return the matrix unitary the implements the qubit clifford C
"""
assert numQubits <= 2, "Oops! I only handle one or two qubits"
if numQubits == 1:
return C1[c]
else:
c = C2Seqs[c]
mat = np.kron(clifford_mat(c[0][0], 1), clifford_mat(c[0][1], 1))
if c[1]:
mat = np.dot(entangling_mat(c[1]), mat)
if c[2]:
mat = np.dot(
np.kron(
clifford_mat(c[2][0], 1), clifford_mat(c[2][1], 1)), mat)
return mat Example 12
def proposals(self, X, return_index=False):
"""Retrieve proposals for a video based on its duration
Parameters
----------
X : ndarray
m x 1 array with video duration
Outputs
-------
Y : ndarray
m * n_prop x 2 array with temporal proposals
"""
if self.priors is None:
raise ValueError('model has not been trained')
Y = np.kron(np.expand_dims(X, 1), self.priors)
idx = np.repeat(np.arange(X.shape[0]), self.n_prop)
if return_index:
return Y, idx
return Y Example 13
def proposals(self, X, return_index=False):
"""Retrieve proposals for a video based on its duration
Parameters
----------
X : ndarray
m x 1 array with video duration
Outputs
-------
Y : ndarray
m * n_prop x 2 array with temporal proposals
"""
if self.priors is None:
raise ValueError('model has not been trained')
Y = np.kron(np.expand_dims(X, 1), self.priors)
idx = np.repeat(np.arange(X.shape[0]), self.n_prop)
if return_index:
return Y, idx
return Y Example 14
def apply_gate(self,gate,on_qubit_name):
on_qubit=self.qubits.get_quantum_register_containing(on_qubit_name)
if len(on_qubit.get_noop()) > 0:
print "NOTE this qubit has been measured previously, there should be no more gates allowed but we are reverting that measurement for consistency with IBM's language"
on_qubit.set_state(on_qubit.get_noop())
on_qubit.set_noop([])
if not on_qubit.is_entangled():
if on_qubit.get_num_qubits()!=1:
raise Exception("This qubit is not marked as entangled but it has an entangled state")
on_qubit.set_state(gate*on_qubit.get_state())
else:
if not on_qubit.get_num_qubits()>1:
raise Exception("This qubit is marked as entangled but it does not have an entangled state")
n_entangled=len(on_qubit.get_entangled())
apply_gate_to_qubit_idx=[qb.name for qb in on_qubit.get_entangled()].index(on_qubit_name)
if apply_gate_to_qubit_idx==0:
entangled_gate=gate
else:
entangled_gate=Gate.eye
for i in range(1,n_entangled):
if apply_gate_to_qubit_idx==i:
entangled_gate=np.kron(entangled_gate,gate)
else:
entangled_gate=np.kron(entangled_gate,Gate.eye)
on_qubit.set_state(entangled_gate*on_qubit.get_state()) Example 15
def cov(M):
"""
Compute the sample covariance matrix of a 2D matrix.
Parameters:
M: `numpy array`
2d matrix of HSI data (N x p)
Returns: `numpy array`
sample covariance matrix
"""
N = M.shape[0]
u = M.mean(axis=0)
M = M - np.kron(np.ones((N, 1)), u)
C = np.dot(M.T, M) / (N-1)
return C Example 16
def com(A, d1, d2):
"""Calculation of the center of mass for spatial components
Inputs:
A: np.ndarray
matrix of spatial components (d x K)
d1: int
number of pixels in x-direction
d2: int
number of pixels in y-direction
Output:
cm: np.ndarray
center of mass for spatial components (K x 2)
"""
nr = np.shape(A)[-1]
Coor = dict()
Coor['x'] = np.kron(np.ones((d2, 1)), np.expand_dims(range(d1), axis=1))
Coor['y'] = np.kron(np.expand_dims(range(d2), axis=1), np.ones((d1, 1)))
cm = np.zeros((nr, 2)) # vector for center of mass
cm[:, 0] = np.dot(Coor['x'].T, A) / A.sum(axis=0)
cm[:, 1] = np.dot(Coor['y'].T, A) / A.sum(axis=0)
return cm Example 17
def find_activity_intervals(C,Npeaks = 5, tB=-5, tA = 25, thres = 0.3):
import peakutils
K,T = np.shape(C)
L = []
for i in range(K):
indexes = peakutils.indexes(C[i,:],thres=thres)
srt_ind = indexes[np.argsort(C[i,indexes])][::-1]
srt_ind = srt_ind[:Npeaks]
L.append(srt_ind)
LOC = []
for i in range(K):
if len(L[i])>0:
interval = np.kron(L[i],np.ones(tA-tB,dtype=int)) + np.kron(np.ones(len(L[i]),dtype=int),np.arange(tB,tA))
interval[interval<0] = 0
interval[interval>T-1] = T-1
LOC.append(np.array(list(set(interval))))
else:
LOC.append(None)
return LOC Example 18
def paulidouble(i, j, tensor=True):
pauli_i = pauli(i)
pauli_j = pauli(j)
outer = np.zeros((4, 4), dtype=np.complex)
outer[:2, :2] = pauli_i
outer[2:, 2:] = pauli_j
if tensor:
outer.shape = (2, 2, 2, 2)
return outer
# def paulitwo_left(i):
# return np.kron(pauli(i), pauli(0))
# def paulitwo_right(i):
# return np.kron(pauli(0), pauli(i))
# def newrho_DK(Lket, theta_ij):
# Lbra = np.conjugate(L_before)
# theta_star = np.conjugate(theta_ij)
# in_bracket = scon(Lbra, Lket, theta_ij, theta_star,
# [1], [1], [2, 3, 1, Example 19
def to_0xyz_basis(ptm):
"""Transform a Pauli transfer in the 0xy1 basis [1],
to the the usual 0xyz. The inverse of to_0xy1_basis.
ptm: The input transfer matrix in 0xy1 basis. Must be 4x4.
[1] Daniel Greenbaum, Introduction to Quantum Gate Set Tomography, http://arxiv.org/abs/1509.02921v1
"""
ptm = np.array(ptm)
if ptm.shape == (4, 4):
trans_mat = basis_transformation_matrix
return np.dot(trans_mat, np.dot(ptm, trans_mat))
elif ptm.shape == (16, 16):
trans_mat = np.kron(basis_transformation_matrix, basis_transformation_matrix)
return np.dot(trans_mat, np.dot(ptm, trans_mat))
else:
raise ValueError("Dimensions wrong, must be one- or two Pauli transfer matrix ") Example 20
def apply_two_ptm(self, bit0, bit1, two_ptm):
"""Apply a two_bit_ptm between bit0 and bit1.
"""
self.ensure_dense(bit0)
self.ensure_dense(bit1)
ptm0 = np.eye(4)
if bit0 in self.single_ptms_to_do:
for ptm2 in self.single_ptms_to_do[bit0]:
ptm0 = ptm2.dot(ptm0)
del self.single_ptms_to_do[bit0]
ptm1 = np.eye(4)
if bit1 in self.single_ptms_to_do:
for ptm2 in self.single_ptms_to_do[bit1]:
ptm1 = ptm2.dot(ptm1)
del self.single_ptms_to_do[bit1]
full_two_ptm = np.dot(two_ptm, np.kron(ptm1, ptm0))
self.full_dm.apply_two_ptm(self.idx_in_full_dm[bit0],
self.idx_in_full_dm[bit1], full_two_ptm) Example 21
def kron_all(op,num,op_2):
# returns an addition of sth like xii + ixi + iix for op =x and op_2 =i
total = np.zeros([len(op)**num,len(op)**num])
a=op
for jj in range(num):
if jj != 0:
a = op_2
else:
a = op
for ii in range(num-1):
if (jj - ii) == 1:
b = op
else:
b = op_2
a = np.kron(a,b)
total = total + a
return a Example 22
def outer_product(input_array):
r'''
Takes a `NumPy` array and returns the outer (dyadic, Kronecker) product
with itself. If `input_array` is a vector :math:`\mathbf{x}`, this returns
:math:`\mathbf{x}\mathbf{x}^T`.
'''
la = len(input_array)
# return outer product as numpy array
return np.kron(input_array, input_array).reshape(la, la)
##############
# Decorators #
##############
# Main decorator for fit functions Example 23
def _upsample_filters(self, filters, rate):
"""Upsamples the filters by a factor of rate along the spatial dimensions.
Args:
filters: [h, w, in_depth, out_depth]. Original filters.
rate: An int, specifying the upsampling rate.
Returns:
filters_up: [h_up, w_up, in_depth, out_depth]. Upsampled filters with
h_up = h + (h - 1) * (rate - 1)
w_up = w + (w - 1) * (rate - 1)
containing (rate - 1) zeros between consecutive filter values along
the filters' spatial dimensions.
"""
if rate == 1:
return filters
# [h, w, in_depth, out_depth] -> [in_depth, out_depth, h, w]
filters_up = np.transpose(filters, [2, 3, 0, 1])
ker = np.zeros([rate, rate])
ker[0, 0] = 1
filters_up = np.kron(filters_up, ker)[:, :, :-(rate-1), :-(rate-1)]
# [in_depth, out_depth, h_up, w_up] -> [h_up, w_up, in_depth, out_depth]
filters_up = np.transpose(filters_up, [2, 3, 0, 1])
self.assertEqual(np.sum(filters), np.sum(filters_up))
return filters_up Example 24
def test_perform(self):
if not imported_scipy:
raise SkipTest('kron tests need the scipy package to be installed')
for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
x = tensor.tensor(dtype='floatX',
broadcastable=(False,) * len(shp0))
a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
if len(shp0) + len(shp1) == 2:
continue
y = tensor.tensor(dtype='floatX',
broadcastable=(False,) * len(shp1))
f = function([x, y], kron(x, y))
b = self.rng.rand(*shp1).astype(config.floatX)
out = f(a, b)
# Newer versions of scipy want 4 dimensions at least,
# so we have to add a dimension to a and flatten the result.
if len(shp0) + len(shp1) == 3:
scipy_val = scipy.linalg.kron(
a[numpy.newaxis, :], b).flatten()
else:
scipy_val = scipy.linalg.kron(a, b)
utt.assert_allclose(out, scipy_val) Example 25
def testCholesky(self):
# Tests the cholesky function
np.random.seed(8)
# generating two symmetric positive-definite tt-cores
L_1 = np.tril(np.random.normal(scale=2., size=(2, 2)))
L_2 = np.tril(np.random.normal(scale=2., size=(3, 3)))
K_1 = L_1.dot(L_1.T)
K_2 = L_2.dot(L_2.T)
K = np.kron(K_1, K_2)
initializer = tensor_train.TensorTrain([K_1[None, :, :, None],
K_2[None, :, :, None]],
tt_ranks=7*[1])
kron_mat = variables.get_variable('kron_mat', initializer=initializer)
init_op = tf.global_variables_initializer()
with self.test_session() as sess:
sess.run(init_op)
desired = np.linalg.cholesky(K)
actual = ops.full(kr.cholesky(kron_mat)).eval()
self.assertAllClose(desired, actual) Example 26
def _iter_contrasts(n_subjects, factor_levels, effect_picks):
""" Aux Function: Setup contrasts """
from scipy.signal import detrend
sc = []
n_factors = len(factor_levels)
# prepare computation of Kronecker products
for n_levels in factor_levels:
# for each factor append
# 1) column vector of length == number of levels,
# 2) square matrix with diagonal == number of levels
# main + interaction effects for contrasts
sc.append([np.ones([n_levels, 1]),
detrend(np.eye(n_levels), type='constant')])
for this_effect in effect_picks:
contrast_idx = _get_contrast_indices(this_effect + 1, n_factors)
c_ = sc[0][contrast_idx[n_factors - 1]]
for i_contrast in range(1, n_factors):
this_contrast = contrast_idx[(n_factors - 1) - i_contrast]
c_ = np.kron(c_, sc[i_contrast][this_contrast])
df1 = matrix_rank(c_)
df2 = df1 * (n_subjects - 1)
yield c_, df1, df2 Example 27
def test_homoskedastic_direct(cov_data, debias):
x, z, eps, sigma = cov_data
cov = HomoskedasticCovariance(x, eps, sigma, sigma, debiased=debias)
k = len(x)
nobs = x[0].shape[0]
big_x = []
for i in range(k):
row = []
for j in range(k):
if i == j:
row.append(x[i])
else:
row.append(np.zeros((nobs, x[j].shape[1])))
big_x.append(np.concatenate(row, 1))
big_x = np.concatenate(big_x, 0)
xeex = big_x.T @ np.kron(sigma, np.eye(nobs)) @ big_x / nobs
xpxi = _xpxi(x)
direct = xpxi @ xeex @ xpxi / nobs
direct = (direct + direct.T) / 2
assert_allclose(np.diag(direct), np.diag(cov.cov))
s = np.sqrt(np.diag(direct))[:, None]
r_direct = direct / (s @ s.T)
s = np.sqrt(np.diag(cov.cov))[:, None]
c_direct = direct / (s @ s.T)
assert_allclose(r_direct, c_direct, atol=1e-5) Example 28
def test_inner_product_short_circuit(data):
y, x, sigma = data
sigma = np.eye(len(x))
efficient = blocked_inner_prod(x, sigma)
nobs = x[0].shape[0]
omega = np.kron(sigma, np.eye(nobs))
k = len(x)
bigx = []
for i in range(k):
row = []
for j in range(k):
if i == j:
row.append(x[i])
else:
row.append(np.zeros((nobs, x[j].shape[1])))
bigx.append(np.hstack(row))
bigx = np.vstack(bigx)
expected = bigx.T @ omega @ bigx
assert_allclose(efficient, expected) Example 29
def test_diag_product(data):
y, x, sigma = data
efficient = blocked_diag_product(x, sigma)
nobs = x[0].shape[0]
omega = np.kron(sigma, np.eye(nobs))
k = len(x)
bigx = []
for i in range(k):
row = []
for j in range(k):
if i == j:
row.append(x[i])
else:
row.append(np.zeros((nobs, x[j].shape[1])))
bigx.append(np.hstack(row))
bigx = np.vstack(bigx)
expected = omega @ bigx
assert_allclose(efficient, expected) Example 30
def test_quantumnumbers_ordinary():
print 'test_quantumnumbers_ordinary'
a=QuantumNumbers('C',([SQN(-1.0),SQN(0.0),SQN(1.0)],[1,1,1]),protocol=QuantumNumbers.COUNTS)
print 'a: %s'%a
print 'copy(a): %s'%copy(a)
print 'deepcopy(a): %s'%deepcopy(a)
b,permutation=QuantumNumbers.kron([a]*2).sort(history=True)
print 'b: ',b
print 'permutation of b:%s'%permutation
print 'b.reorder(permutation,protocol="EXPANSION"): \n%s'%b.reorder(permutation,protocol="EXPANSION")
print 'b.reorder([4,3,2,1,0],protocol="CONTENTS"): \n%s'%b.reorder([4,3,2,1,0],protocol="CONTENTS")
c=b.to_ordereddict(protocol=QuantumNumbers.COUNTS)
print 'c(b.to_ordereddict(protocol=QuantumNumbers.COUNTS)):\n%s'%('\n'.join('%s: %s'%(key,value) for key,value in c.iteritems()))
print 'QuantumNumbers.from_ordereddict(SQN,c,protocol=QuantumNumbers.COUNTS):\n%s'%QuantumNumbers.from_ordereddict(SQN,c,protocol=QuantumNumbers.COUNTS)
d=b.to_ordereddict(protocol=QuantumNumbers.INDPTR)
print 'd(b.to_ordereddict(protocol=QuantumNumbers.INDPTR)):\n%s'%('\n'.join('%s: %s'%(key,value)for key,value in d.iteritems()))
print 'QuantumNumbers.from_ordereddict(SQN,d,protocol=QuantumNumbers.INDPTR):\n%s'%QuantumNumbers.from_ordereddict(SQN,d,protocol=QuantumNumbers.INDPTR)
print Example 31
def hmm_trans_matrix(self):
# NOTE: more general version, allows different delays, o/w we could
# construct with np.kron
if self._hmm_trans_matrix is None:
ps, delays = map(np.array,zip(*[(d.p,d.delay) for d in self.dur_distns]))
starts, ends = cumsum(delays,strict=True), cumsum(delays,strict=False)
trans_matrix = self._hmm_trans_matrix = np.zeros((ends[-1],ends[-1]))
for (i,j), Aij in np.ndenumerate(self.trans_matrix):
block = trans_matrix[starts[i]:ends[i],starts[j]:ends[j]]
if i == j:
block[:-1,1:] = np.eye(block.shape[0]-1)
block[-1,-1] = 1-ps[i]
else:
block[-1,0] = ps[j]*Aij
return self._hmm_trans_matrix Example 32
def _initR(X, A, lmbdaR):
_log.debug('Initializing R (SVD) lambda R: %s' % str(lmbdaR))
rank = A.shape[1]
U, S, Vt = svd(A, full_matrices=False)
Shat = kron(S, S)
Shat = (Shat / (Shat ** 2 + lmbdaR)).reshape(rank, rank)
R = []
ep = 1e-9
for i in range(len(X)): # parallelize
Rn = Shat * dot(U.T, X[i].dot(U))
Rn = dot(Vt.T, dot(Rn, Vt))
negativeVal = Rn.min()
Rn.__iadd__(-negativeVal+ep)
# if Rn.min() < 0 :
# print("Negative Rn!")
# raw_input("Press Enter: ")
# Rn = np.eye(rank)
# Rn = dot(A.T,A)
R.append(Rn)
print('Initialized R')
return R
# ------------------ Update V ------------------------------------------------ Example 33
def concurrence(self, rho): """Compute the concurrence of the density matrix. :param numpy_array rho: Density matrix :return: The concurrence, see https://en.wikipedia.org/wiki/Concurrence_(quantum_computing). :rtype: complex """ rhoTilde = np.dot( np.dot(np.kron(self.PAULI[1], self.PAULI[1]) , rho.conj()) , np.kron(self.PAULI[1], self.PAULI[1])) rhoRoot = scipy.linalg.fractional_matrix_power(rho, 1/2.0) R = scipy.linalg.fractional_matrix_power(np.dot(np.dot(rhoRoot,rhoTilde),rhoRoot),1/2.0) eigValues, eigVecors = np.linalg.eig(R) sortedEigValues = np.sort(eigValues) con = sortedEigValues[3]-sortedEigValues[2]-sortedEigValues[1]-sortedEigValues[0] return np.max([con,0])
Example 34
def pauli_mpp(nr_sites, local_dim):
r"""Pauli POVM tensor product as MP-POVM
The resulting MP-POVM will contain all tensor products of the
elements of the local Pauli POVM from :func:`mpp.pauli_povm()`.
:param int nr_sites: Number of sites of the returned MP-POVM
:param int local_dim: Local dimension
:rtype: MPPovm
For example, for two qubits the `(1, 3)` measurement outcome is
`minus X` on the first and `minus Y` on the second qubit:
>>> nr_sites = 2
>>> local_dim = 2
>>> pauli = pauli_mpp(nr_sites, local_dim)
>>> xy = np.kron([1, -1], [1, -1j]) / 2
>>> xyproj = np.outer(xy, xy.conj())
>>> proj = pauli.get([1, 3], astype=mp.MPArray) \
... .to_array_global().reshape((4, 4))
>>> abs(proj - xyproj / 3**nr_sites).max() <= 1e-10
True
The prefactor `1 / 3**nr_sites` arises because X, Y and Z are in a
single POVM.
"""
from mpnum.povm import pauli_povm
return MPPovm.from_local_povm(pauli_povm(local_dim), nr_sites) Example 35
def mkron(*args):
"""np.kron() with an arbitrary number of n >= 1 arguments"""
if len(args) == 1:
return args[0]
return mkron(np.kron(args[0], args[1]), *args[2:]) Example 36
def test_partialdot(nr_sites, local_dim, rank, rgen, dtype):
# Only for at least two sites, we can apply an operator to a part
# of a chain.
if nr_sites < 2:
return
part_sites = nr_sites // 2
start_at = min(2, nr_sites // 2)
mpo = factory.random_mpa(nr_sites, (local_dim, local_dim), rank,
randstate=rgen, dtype=dtype)
op = mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
mpo_part = factory.random_mpa(part_sites, (local_dim, local_dim), rank,
randstate=rgen, dtype=dtype)
op_part = mpo_part.to_array_global().reshape((local_dim**part_sites,) * 2)
op_part_embedded = np.kron(
np.kron(np.eye(local_dim**start_at), op_part),
np.eye(local_dim**(nr_sites - part_sites - start_at)))
prod1 = np.dot(op, op_part_embedded)
prod2 = np.dot(op_part_embedded, op)
prod1_mpo = mp.partialdot(mpo, mpo_part, start_at=start_at)
prod2_mpo = mp.partialdot(mpo_part, mpo, start_at=start_at)
prod1_mpo = prod1_mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
prod2_mpo = prod2_mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
assert_array_almost_equal(prod1, prod1_mpo)
assert_array_almost_equal(prod2, prod2_mpo)
assert prod1_mpo.dtype == dtype
assert prod2_mpo.dtype == dtype Example 37
def state_from_string(qubit_state_string):
if not all(x in '01' for x in qubit_state_string):
raise Exception("Description must be a string in binary")
state=None
for qubit in qubit_state_string:
if qubit=='0':
new_contrib=State.zero_state
elif qubit=='1':
new_contrib=State.one_state
if state==None:
state=new_contrib
else:
state=np.kron(state,new_contrib)
return state Example 38
def separate_state(qubit_state):
"""This only works if the state is fully separable at present
Throws exception if not a separable state"""
n_entangled=QuantumRegister.num_qubits(qubit_state)
if list(qubit_state.flat).count(1)==1:
separated_state=[]
idx_state=list(qubit_state.flat).index(1)
add_factor=0
size=qubit_state.shape[0]
while(len(separated_state)<n_entangled):
size=size/2
if idx_state<(add_factor+size):
separated_state+=[State.zero_state]
add_factor+=0
else:
separated_state+=[State.one_state]
add_factor+=size
return separated_state
else:
# Try a few naive separations before giving up
cardinal_states=[State.zero_state,State.one_state,State.plus_state,State.minus_state,State.plusi_state,State.minusi_state]
for separated_state in itertools.product(cardinal_states, repeat=n_entangled):
candidate_state=reduce(lambda x,y:np.kron(x,y),separated_state)
if np.allclose(candidate_state,qubit_state):
return separated_state
# TODO: more general separation methods
raise StateNotSeparableException("TODO: Entangled qubits not represented yet in quantum computer implementation. Can currently do manual calculations; see TestBellState for Examples") Example 39
def apply_two_qubit_gate_CNOT(self,first_qubit_name,second_qubit_name):
""" Should work for all combination of qubits"""
first_qubit=self.qubits.get_quantum_register_containing(first_qubit_name)
second_qubit=self.qubits.get_quantum_register_containing(second_qubit_name)
if len(first_qubit.get_noop())>0 or len(second_qubit.get_noop())>0:
raise Exception("Control or target qubit has been measured previously, no more gates allowed")
if not first_qubit.is_entangled() and not second_qubit.is_entangled():
combined_state=np.kron(first_qubit.get_state(),second_qubit.get_state())
if first_qubit.get_num_qubits()!=1 or second_qubit.get_num_qubits()!=1:
raise Exception("Both qubits are marked as not entangled but one or the other has an entangled state")
new_state=Gate.CNOT2_01*combined_state
if State.is_fully_separable(new_state):
second_qubit.set_state(State.get_second_qubit(new_state))
else:
self.qubits.entangle_quantum_registers(first_qubit,second_qubit)
first_qubit.set_state(new_state)
else:
if not first_qubit.is_entangled_with(second_qubit):
# Entangle the state
combined_state=np.kron(first_qubit.get_state(),second_qubit.get_state())
self.qubits.entangle_quantum_registers(first_qubit,second_qubit)
else:
# We are ready to do the operation
combined_state=first_qubit.get_state()
# Time for more meta programming!
# Select gate based on indices
control_qubit_idx,target_qubit_idx=first_qubit.get_indices(second_qubit)
gate_size=QuantumRegister.num_qubits(combined_state)
try:
exec 'gate=Gate.CNOT%d_%d%d' %(gate_size,control_qubit_idx,target_qubit_idx)
except:
print 'gate=Gate.CNOT%d_%d%d' %(gate_size,control_qubit_idx,target_qubit_idx)
raise Exception("Unrecognized combination of number of qubits")
first_qubit.set_state(gate*combined_state) Example 40
def computeMisfit2(self,q):
assert self.r0 is not None, "Must have current estimate of polarizations"
assert self.dunc is not None, "Must have set uncertainties"
assert self.dobs is not None, "Must have observed data"
dunc = mkvc(self.dunc)
dobs = mkvc(self.dobs)
r0 = self.r0
Hp = self.computeHp(r0=r0,update=False)
Brx = self.computeBrx(r0=r0,update=False)
P = self.computeP(Hp,Brx)
N = np.size(dobs)
M = len(self.times)
Px = np.kron(np.diag(np.ones(M)),P)
dpre = np.dot(Px,q)
v = (dpre-dobs)/dunc
Phi = np.dot(v.T,v)
return Phi/N Example 41
def updatePolarizationsQP(self,r0,UB):
# Set operator and solution array
Hp = self.computeHp(r0=r0)
Brx = self.computeBrx(r0=r0)
P = self.computeP(Hp,Brx)
dunc = self.dunc
dobs = self.dobs
K = np.shape(dobs)[1]
q = np.zeros((6,K))
# Bounds
ub = UB*np.ones(6)
e = matrix(np.r_[np.zeros(9),ub])
# Constraints
C1 = np.r_[np.c_[-1,0,0,0,0,0],np.c_[0,0,0,-1,0,0],np.c_[0,0,0,0,0,-1]]
C2 = np.r_[np.c_[-0.5,1,0,-0.5,0,0],np.c_[-0.5,0,1,0,0,-0.5],np.c_[0,0,0,-0.5,1,-0.5]]
C3 = np.r_[np.c_[-0.5,-1,0,-0.5,0,0],np.c_[-0.5,0,-1,0,0,-0.5],np.c_[0,0,0,-0.5,-1,-0.5]]
# G = sp.kron(sp.diags([-np.ones(K),np.ones(K)],[0,1],shape=(K-1,K)),sp.diags(np.ones(6)))
I = np.diag(np.ones(6))
G = np.r_[C1,C2,C3,I]
G = matrix(G)
solvers.options['show_progress'] = False
for kk in range(0,K):
LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
RHS = dobs[:,kk]/dunc[:,kk]
A = matrix(np.dot(LHS.T,LHS))
b = matrix(np.dot(LHS.T,-RHS))
sol = solvers.qp(A, b, G=G, h=e)
q[:,kk] = np.reshape(sol['x'],(6))
self.q = q Example 42
def updatePolarizationsQP(self,r0,UB):
# Set operator and solution array
Hp = self.computeHp(r0=r0)
Brx = self.computeBrx(r0=r0)
P = self.computeP(Hp,Brx)
dunc = self.dunc
dobs = self.dobs
K = np.shape(dobs)[1]
q = np.zeros((6,K))
# Bounds
ub = UB*np.ones(6)
e = matrix(np.r_[np.zeros(9),ub])
# Constraints
C1 = np.r_[np.c_[-1,0,0,0,0,0],np.c_[0,0,0,-1,0,0],np.c_[0,0,0,0,0,-1]]
C2 = np.r_[np.c_[-0.5,1,0,-0.5,0,0],np.c_[-0.5,0,1,0,0,-0.5],np.c_[0,0,0,-0.5,1,-0.5]]
C3 = np.r_[np.c_[-0.5,-1,0,-0.5,0,0],np.c_[-0.5,0,-1,0,0,-0.5],np.c_[0,0,0,-0.5,-1,-0.5]]
# G = sp.kron(sp.diags([-np.ones(K),np.ones(K)],[0,1],shape=(K-1,K)),sp.diags(np.ones(6)))
I = np.diag(np.ones(6))
G = np.r_[C1,C2,C3,I]
G = matrix(G)
solvers.options['show_progress'] = False
for kk in range(0,K):
LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
RHS = dobs[:,kk]/dunc[:,kk]
A = matrix(np.dot(LHS.T,LHS))
b = matrix(np.dot(LHS.T,-RHS))
sol = solvers.qp(A, b, G=G, h=e)
q[:,kk] = np.reshape(sol['x'],(6))
self.q = q Example 43
def area(self):
"""Face areas"""
if getattr(self, '_area', None) is None:
if self.nCy > 1:
raise NotImplementedError('area not yet implemented for 3D '
'cyl mesh')
areaR = np.kron(self.hz, 2*pi*self.vectorNx)
areaZ = np.kron(np.ones_like(self.vectorNz), pi*(self.vectorNx**2 -
np.r_[0, self.vectorNx[:-1]]**2))
self._area = np.r_[areaR, areaZ]
return self._area Example 44
def vol(self):
"""Volume of each cell"""
if getattr(self, '_vol', None) is None:
if self.nCy > 1:
raise NotImplementedError('vol not yet implemented for 3D '
'cyl mesh')
az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2)
self._vol = np.kron(self.hz, az)
return self._vol
####################################################
# Operators
#################################################### Example 45
def aveF2CC(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CC', None) is None:
n = self.vnC
if self.isSymmetric:
avR = utils.av(n[0])[:, 1:]
avR[0, 0] = 1.
self._aveF2CC = ((0.5)*sp.hstack((sp.kron(utils.speye(n[2]),
avR),
sp.kron(utils.av(n[2]),
utils.speye(n[0]))),
format="csr"))
else:
raise NotImplementedError('wrapping in the averaging is not '
'yet implemented')
# self._aveF2CC = (1./3.)*sp.hstack((utils.kron3(utils.speye(n[2]),
# utils.speye(n[1]),
# utils.av(n[0])),
# utils.kron3(utils.speye(n[2]),
# utils.av(n[1]),
# utils.speye(n[0])),
# utils.kron3(utils.av(n[2]),
# utils.speye(n[1]),
# utils.speye(n[0]))),
# format="csr")
return self._aveF2CC Example 46
def normalize2D(x):
return x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2)) Example 47
def normalize3D(x):
return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) Example 48
def test_kron_matrix(self, level=rlevel):
# Ticket #71
x = np.matrix('[1 0; 1 0]')
assert_equal(type(np.kron(x, x)), type(x)) Example 49
def kron(a, b):
"""Returns the kronecker product of two arrays.
Args:
a (~cupy.ndarray): The first argument.
b (~cupy.ndarray): The second argument.
Returns:
~cupy.ndarray: Output array.
.. seealso:: :func:`numpy.kron`
"""
a_ndim = a.ndim
b_ndim = b.ndim
if a_ndim == 0 or b_ndim == 0:
return cupy.multiply(a, b)
ndim = b_ndim
a_shape = a.shape
b_shape = b.shape
if a_ndim != b_ndim:
if b_ndim > a_ndim:
a_shape = (1,) * (b_ndim - a_ndim) + a_shape
else:
b_shape = (1,) * (a_ndim - b_ndim) + b_shape
ndim = a_ndim
axis = ndim - 1
out = core.tensordot_core(a, b, None, a.size, b.size, 1, a_shape + b_shape)
for _ in six.moves.range(ndim):
out = core.concatenate_method(out, axis=axis)
return out Example 50
def selftest1():
# Generate data
D0 = 5
K1 = 2
K2 = 2
NperClass = 500
N = NperClass*K1*K2
#l = 1.0e-3
X = np.zeros((D0,NperClass,K1,K2))
y1 = np.zeros((NperClass,K1,K2),dtype=int)
y2 = np.zeros((NperClass,K1,K2),dtype=int)
bias1 = np.random.normal(scale=1.0,size=(D0,K1))
bias2 = np.random.normal(scale=1.0,size=(D0,K2))
for k1 in range(K1):
for k2 in range(K2):
X[:,:,k1,k2] = \
np.random.normal(scale=0.25, size=(D0,NperClass)) \
+ np.kron(np.ones((1,NperClass)),bias1[:,k1].reshape((D0,1))) \
+ np.kron(np.ones((1,NperClass)),bias2[:,k2].reshape((D0,1)))
y1[:,k1,k2] = k1*np.ones((NperClass,))
y2[:,k1,k2] = k2*np.ones((NperClass,))
X = X.reshape((D0,N))
y1 = y1.reshape((N,))
y2 = y2.reshape((N,))
E,dd = run(X,y1,y2)
print dd
plt.figure(1)
plt.clf()
plt.subplot(1,2,1)
plt.plot(dd)
plt.subplot(1,2,2)
plt.imshow(E, aspect='auto', interpolation='none')
plt.colorbar()
plt.show() 